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The  low-level  interpretation  of  images- provides  constraints  on  3D  surface  shape  at  multiple  reso¬ 
lutions,  but  typically  only  at  scattered  locations  over  the  visual  field.  Subsequent  visual  processing 
can  be  facilitated  substantially  if  die  scattered  shape  constraints  are  immediately  transformed  into 
visible-surface  representations  that  unambiguously  specify  surface  shape  at  every  image  point.  The 
required  transformation  is  shown  to  lead  to  an  ill-posed  surface  reconstruction  problem.  A  well- 
posed  variational  principle  formulation  is  obtained  by  invoking  “controlled  continuity,"  a  physically 
nomostrictivc  (generic)  assumption  about  surfaces  which  is  nonetheless  strong  enough  to  guarantee 
unique  solutions.  The  variational  principle,  which  admits  an  appealing  physical  interpretation,  is 
locally  discretized  by  applying  the  finite  element  method  to  a  piecewise,  finite  element  represen¬ 
tation  of  surfaces.  This  forms  tire  mathematical  basis  of  a  unified  and  general  framework  for 
computing  visible-surface  representations.  The  computational  framework  unifies  formal  solutions 
to  the  key  problems  of  (i)  integrating  multiscale  constraints  on  surface  depth  and  orientation  from 
multiple  visual  sources,  (ii)  interpolating  these  scattered  constraints  into  dense,  piecewise  smooth 
surfaces,  (iii)  discovering  surface  depth  and  orientation  discontinuities  and  allowing  them  to  restrict 
interpolation  appropriately,  and  (iv)  overcoming  the  immense  computational  burden  of  fine  resolu¬ 
tion  surface  reconstruction.  An  efficient  surface  reconstruction  algorithm  is  developed.  If  exploits 
multircsolution  hierarchies  of  cooperative  relaxation  processes  and  is  suitable  for  implementation  on 
massively  parallel  networks  of  simple,  locally  interconnected  processors.  The  algorithm  is  evaluated 
empirically  in  a  diversity  of  applications. 
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formulation  is  obtained  by  invoking  "controlled  continuity,"  a  physically  non-restictive 
(generic)  assumption  about  surfaces  which  is  nonetheless  strong  enough  to  guarantee 
unique  solutions.  The  variational  principle,  which  admits  an  appealing  physical 
interpretation,  is  locally  discretized  by  applying  the  finite  element  method 
to  a  piecewise,  finite  element  representation  of  surfaces.  This  forms  the  mathematical 
basis  of  a  unified  and  general  framework  for  computingvisible-surface  representations. 
The  computational  framework  unifies  formal  solutions  to  the  key  problems  of 
(i)  intergrating  multiscale  constraints  on  surface  depth  and  orientation  from 
multiple  visual  sources,  (ii)  interpolating  these  scattered  constraints  into 
dense,  piecewise  smooth  surfaces,  (iii)  discovering  surface  depth  and  orientation 
discontinuities  and  allowing  them  to  restrict  interpolation  appropiately  ,  and 
(iv)  overcoming  the  immense  computational  burden  of  fine  resolution  surface 
reconstruction.  An  efficient  surface  reconstruction  algorithm  is  developed. 

It  exploits  multiresolution  hierarchies  of  cooperative  relaxation  processes 
and  is  suitable  for  implementation  on  massively  parallel  networks  of  simple, 
locally  interconnected  processors.  The  algorithm  is  evaluated  empirically  in 
a  diversity  of  applications. 
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Computing  Visible-Surface  Representations 


Terzopoulos 

1.  Introduction 


Over  thirty  years  ago,  J.J.  Gibson  [1950J  made  the  seminal  conjecture  that  natural  human  peteption 
amounts  to  the  perception  of  visible  surfaces.  The  explicit  representation  of  visible  surfaces,  an 
intermediate  goal  of  computational  vision,  has  since  attracted  considerable  interest. 

The  computational  framework  offered  in  this  paper  addresses,  in  a  unified  way,  certain  visual 
information  processing  tasks  involved  in  the  representation  of  visible  surfaces.  Particular  emphasis 
is  placed  on  utilizing  highly  parallel,  cooperative  processing  to  integrate  surface  shape  information 
over  multiple  visual  sources,  to  fuse  it  across  a  multiplicity  of  spatial  resolutions,  and  to  maintain  the 
global  consistency  of  the  resulting  distributed  shape  representations.  The  issues  are  first  investigated 
in  terms  of  a  surface  reconstruction  model  rooted  in  mathematical  physics.  This  formal  analysis  is 
augmemed  by  an  empirical  study  of  the  resulting  algorithms,  which  feature  multircsolution  iterative 
processing  within  hierarchical  surface  shape  representations.  The  approach  is  guided  by  current 
knowlcc  gc  of  how  humans  perceive  visible  surfaces,  while  applications  in  machine  vision  provide 
a  testbed  for  the  algorithms,  j  LjZ 

The  remainder  of  this  introdfactory  section  examines  the  role  of  surface  representations  in  early  ) 

visual  processing,  outlines  tire  key  computational  problems  that  will  be  of  primary  concern,  and 
reviews  some  relevant  prior  work.  ^  - * - — 


1.1.  Ecrly  Visual  Processing  and  Visible-Surface  Representations 


Early  vision  comprises  a  set  of  processes  which  specialize  in  recovering  the  physical  properties  of 
visible  surfaces  in  a  3D  scene  from  2D  images  of  the  scene.  rrhcy  apply  generic  assumptions  about 
the  physical  world  and  the  imaging  process  to  infer  3D  surface  shape  constraints  by  interpreting 
specific  image  cues,  such  as  stereoscopic  disparity,  motion,  texture,  contours,  and  shading.  These 
conceptually  independent  shape  estimation  processes  fall  into  two  broad  categories. 

The  first  category  comprises  what  arc  commonly  referred  to  as  correspondence  processes.  They 
operate  over  multiple  image  frames  of  a  scene  taken  across  space  or  over  time.  Paradigm  examples 
are  stercopsis  and  structure  from  motion  (sec,  e.g.,  the  review  articles  [Poggio  and  Poggio,  1984] 
and  [Ullman,  1983]).  Stercopsis  is  driven  by  computations  on  typically  two  image  frames  taken 
simultaneously,  but  from  different  spatial  positions.  The  basic  structure  from  motion  computation 
involves  frames  taken  from  the  same  position,  but  at  different  times.  If  correspondences  can  be 
established  across  the  frames,  between  image  features  which  originate  from  the  same  point  on 
a  visible  surface  (not  a  trivial  problem),  then  the  depth  (i.e.,  3D  distance)  to  such  points  can 
be  estimated  by  triangulation,  given  the  disparity  (i.e.,  2D  displacement)  between  corresponding 
features  as  well  as  some  knowledge  of  the  imaging  geometry. 

The  second  category  of  shape  estimation  processes  involve  computations  on  a  single  static 
frame.  Perspective  projection  of  3D  scenes  onto  images  imparts  a  systematic  distortion  to  imaged 
surface  properties  such  as  shading,  texture,  and  contours.  A  major  part  of  this  distortion  can  be 
attributed  to  the  relative  orientations  of  visible  surfaces  with  respect  to  the  viewer.  In  principle, 
it  is  possible  to  estimate  surface  orientation  by  measuring  and  interpreting  such  distortion  in  the 
image,  litis  is  the  basis  of  practical  approaches  to  recovering  surface  shape  from  shading,  texture, 
and  contours  [Ikcuchi  and  Horn,  1981;  Horn  and  Brooks,  1985;  Kcndcr,  1980;  Witkin,  1981;  Brady 
and  Yuillc,  1984]. 

The  combined  output  of  the  shape  estimation  processes  is  best  collected  into  intennediate 
representations  of  the  3D  shapes  and  configurations  of  visible  surfaces,  which  we  will  refer  to 
as  visible-surface  representations.  Notable  among  proposed  visible-surface  representations  arc  the 
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depth  and  needle  maps  of  Horn  [1982],  the  intrinsic  images  of  Barrow  and  Tcncnbaum  [1978],  and 
the  2^-1)  sketch  of  Marr  and  Nishihara  [1978].  For  humans,  the  perception  of  visible  surfaces  is 
generally  immediate,  involuntary,  and  seems  to  precede  (object)  recognition,  'litis  strongly  suggests 
the  existence  of  a  visual  process  that  autonomously  computes  visible-surface  representations.  Aside 
from  the  perceptual  evidence,  the  availability  of  explicit  visible-surface  representations  can  also 
substantially  facilitate  subsequent  surface  analysis  tasks  in  machine  vision. 

Since  early  visual  processing  provides  relative  surface  shape  estimates  with  respect  to  the 
viewer,  it  is  most  natural  to  define  the  basic  shape  primitives  of  visible-surface  representations  in  a 
viewer-centered  coordinate  system.  Moreover,  the  primitives  should  be  computationally  compatible 
with  the  local  depth  and  orientation  measurements  (as  well  as  discontinuities)  that  arc  provided  by 
the  various  shape  estimation  processes.  These  criteria  arc  satisfied  by  a  particularly  appealing  class 
of  local,  piecewise  shape  primitives  known  as  finite  elements. 

A  crucial  realization  is  that  shape  estimates  can  be  provided  at  multiple  resolutions.  Indeed, 
multircsolution  spatial  frequency  channels  have  been  identified  psychophysically  in  the  human 
visual  system  (c.g.,  [Braddick  el  al.,  1978]).  Their  existence  has  influenced  the  design  of  early 
visual  algorithms  (e.g„  [Marr,  1982]).  In  addition,  machine  vision  research  has  demonstrated  that 
multircsolution  processing  effectively  bridges  fine  and  coarse  image  structure,  while  it  simultaneously 
increases  computational  efficiency  (e.g.,  [Roscnfcld.  1984]).  Hence,  a  multiresolution  organization 
of  visible-surface  representations  is  most  desirable  [Tcrzopoulos,  1982,  1983a], 

1.2.  Key  Problems  of  Visible-Surface  Reconstruction 

The  main  topic  of  concern  in  this  paper  is  the  development  of  a  visible-surface  reconstruction  process 
responsible  for  generating  and  dynamically  maintaining  visible  surface  representations.  Whether 
the  intention  is  to  model  human  vision  or  to  design  competent  artificial  vision  systems,  this  process 
must  solve  four  key  problems  [Tcrzopoulos,  1983b.  1984]:  (i)  the  constraint  integration  problem, 
(ii)  the  interpolation  problem,  (iii)  the  discontinuity  problem,  and  (iv)  the  computational  efficiency 
problem.  We  elaborate  on  each  of  these  problems  next. 

(i)  The  Constraint  Integration  Problem :  Each  specialized  visual  process  may  be  thought  of  as  a 
quasi-independent  source  of  information  partially  constraining  the  shapes  of  visible  surfaces. 
The  human  visual  system  is  reliable  and  robust  because  it  integrates  the  various  processes, 
enabling  them  to  complement  one  another.  The  integration  of  multiple  sources  of  information 
introduces  redundancy,  which  is  necessary  not  only  to  resolve  potential  ambiguities,  but  also 
to  overcome  the  detrimental  effects  of  noise  and  inaccuracies  in  the  initial  shape  estimates. 
T  he  constraint  integration  problem  is  fundamentally  one  of  devising  an  effective  means  of 
integrating  all  available  surface  depth  and  orientation  constraints  (and  discontinuities)  within  a 
cooperative  visible-surface  reconstruction  process. 

(ii)  The  Interpolation  Problem:  It  is  widely  accepted  that  initial  descriptions  of  images  ought  to 
make  explicit  the  occurrence  and  local  2D  structure  of  image  features  that  arc  correlated  to 
salient  events  on  physical  surfaces  (markings,  boundaries,  etc.).  This  is  the  essence,  for  instance, 
of  Marr’s  “primal  sketch”  representation  of  significant  image  irradiancc  changes  (edges)  [Marr, 
1982].  Generally,  such  salient  features  do  not  occur  everywhere  over  the  visual  field.  The 
initial  representation  of  images  as  a  sparse  set  of  features  implies  that  surface  shape  constraints 
generated  by  the  specialized  processes  will  also  be  scattered  over  a  subset  of  image  points.  It 
is  fascinating,  however,  that  the  human  visual  system  systematically  interprets  visual  stimuli 
such  as  sparse  random  dot  stereograms  as  coherent  31)  surfaces  [Jules/,  1971],  Indeed,  these 
stereograms  continue  to  elicit  perceptions  of  dense  surfaces,  even  when  die  density  of  dots 
carrying  disparity  information  has  been  reduced  until  depth  is  unspecified  over  98  percent  of 
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the  visible  surface  area  (sec  Fig.  1).  It  therefore  appears  that  the  surface  reconstruction  process 
is  smoothly  ‘'filling  in  the  gaps.”  This  phenomenon  has  been  the  subject  of  some  psychophysical 
investigation  (c.g..  {Collett,  1984]).  The  interpolation  problem  of  visible-surface  rcconst notion 
challenges  us  to  devise  a  scheme,  consistent  with  human  perception,  for  propagating  shape 
information  into  indeterminate  regions  (devoid  of  shape  estimates)  from  places  where  it  is 
available. 


Figure  J.  A  sparse  random  dot  stereogram.  Binocular  fusion  of  this  stereogram  reveals  a  planar  surface  as  a 
central,  opaque,  textured  square  suspended  nearer  in  depth  over  a  similarly  textured  background.  Vivid  depth 
discontinuities  separate  the  dense  surfaces. 


(iii)  The  Discontinuity  Problem:  Visual  discontinuities  result  from  significant,  spatially-localized 
changes  in  the  physical  world,  particularly  abrupt  changes  in  surface  structure.  Both  depth 
and  orientation  discontinuities  arc  perceptually  relevant  and  provide  vital  boundary  conditions 
for  surface  reconstruction.  Discontinuities  in  depth  occur  at  occluding  contours,  along  which 
a  surface  in  the  scene  occludes  itself  or  another  surface.  Orientation  discontinuities  occur  at 
creases  or  cusps  of  an  otherwise  continuous  surface.  In  addition  to  the  perception  of  coherent 
surfaces,  random  dot  stereograms  elicit  vivid  perceptions  of  surface  discontinuities  at  abrupt 
disparity  changes  (sec  Fig.  1).  The  discontinuity  problem  amounts  to  (1)  finding  both  depth 
and  orientation  discontinuities  in  surfaces,  and  (2)  dealing  with  their  presence  during  visible- 
surface  reconstruction;  i.c.,  allowing  discontinuities  to  limit  the  otherwise  smooth  interpolation 
of  shape  constraints. 

(iv)  The  Computational  Efficiency  Problcnu  Visible-surface  reconstruction  at  the  resolution  of  the 
image  imposes  an  immense  computational  burden  on  both  biological  and  artificial  vision  sys¬ 
tem..  Nevertheless,  visible-surface  representations  must  be  computed  quickly  if  they  arc  to 
be  of  any  practical  value.  It  is  generally  accepted  that  to  achieve  the  necessary  performance, 
visit  il  algorithms  and  mechanisms  must  cinphasi/e  parallelism  (Ballard  et  a /.,  1983];  however, 
visible-surface  reconstruction  is  compute  bound  to  the  point  where  the  fundamental  limitations 
of  massively  parallel  mechanisms,  particularly  with  respect  to  global  inlcrproccssor  commu¬ 
nications,  lead  to  severe  inellicicncies.  I  he  computational  efficiency  problem  is  to  develop  a 
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visible-surface  reconstruction  process  that  not  only  exploits  parallelism,  but  also  overcomes  co¬ 
operative  communication  bottlenecks  to  compute  visible-surface  representations  quickly,  given 
suitable  architectures.  Kor  the  reasons  outlined  in  the  foregoing  section  and  following  our 
previous  work  (Terzopoulos,  1982,  1983a],  our  solution  to  this  problem  hinges  on  the  idea  of 
multircsolution  structuring  of  visual  representations  and  associated  cooperative  processes. 

1.3.  Prior  Work 

There  has  been  some  prior  work  relating  to  the  computation  of  visible-surface  representations. 
Barrow  and  Tcnenbaum  [1979]  describe  an  approach  to  reconstructing  smooth  surfaces  from  noisy 
visual  data.  1'his  approach  did  not  apply  to  general  classes  of  surfaces,  however,  and  the  proposed 
relaxation  algorithms  were  not  supported  by  a  firm  mathematical  analysis.  Nevertheless,  Barrow 
and  Tcncnbaum’s  (1978]  basic  model  of  intrinsic  images  and  much  of  the  philosophy  underlying 
their  computation  seems  appropriate,  and  it  has  influenced  our  approach. 

The  interpolation  problem  is  related  to  classical  spline  approximation.  A  number  of  well-known 
surface  approximation  methods  for  scattered  data  are  reviewed  by  Schumakcr  (1976],  Grimson 
(1983]  employed  one  of  these  methods  for  the  continuous  interpolation  of  visual  surfaces  from  depth 
constraints;  a  minimization  scheme  involving  a  particular  functional  containing  second  derivatives 
(he  referred  to  it  as  the  “quadratic  variation”).  Brady  and'  Horn  [1983]  observe  that  this  functional 
is  related  to  the  bending  energy  of  a  thin  plate  (a  connection  noted  by  Duchon  [1977]),  and  the 
thin  plate  model  was  developed  further  by  Terzopoulos  [1983a]  (see  also  (Blake,  1984]). 

Interestingly,  thin  plate  intcrpolants  have  appeared  in  other  areas,  including  the  interpolation 
of  aircraft  wing  deflections  [Harder  and  Dcsmarais,  1972],  interpolation  of  meteorological  fields 
[Wahba  and  Wcndclbergcr.  1980],  and  the  interpolation  of  digital  terrain  maps  [Briggs,  1974; 
Bolondi  el  ai,  1976],  In  this  latter  paper  there  is  some  concern  for  the  presence  of  discontinuities 
(faults). 

Following  Ullman  [1979]  and  others,  Grimson  [1983]  pursued  "biologically  feasible,”  parallel 
and  iterative  algorithms  for  surface  interpolation.  A  serious  drawback  of  algorithms  which  satisfy 
these  criteria  is  that  they  often  converge  excruciatingly  slowly  for  problems  of  reasonable  size. 
The  idea  of  multircsolution  surface  reconstruction  exploiting  multigrid  relaxation  methods  was 
shown  to  overcome  this  problem  while  adhering  to  biological  feasibility  (Terzopoulos,  1982,  1983a]. 
Ihc  multiresolution  methodology  yields  efficient  algorithms  not  only  for  the  surface  reconstruction 
problem  but  for  other  visual  problems  as  well  [Terzopoulos,  1984]. 

In  retrospect,  although  progress  has  been  made,  a  satisfactory  computational  theory  of  visible- 
surface  representations  has  been  elusive.  This  is  largely  a  consequence  of  the  significant  technical 
obstacles  encountered  in  devising  formal  solutions  to  all  four  key  problems  of  visible-surface 
reconstruction  within  a  unified  computational  framework.  Ihc  difficulty  of  the  task  appears  to  have 
evoked  some  skepticism  as  to  the  actual  computability  (hence,  even  the  usefulness)  of  intrinsic 
surface  representations  [Witkin  and  Tcnenbaum,  1983],  Based  on  the  theoretical  generality  of 
our  approach  and  the  accompanying  empirical  results,  however,  we  believe  such  skepticism  to  be 
premature. 


2.  Mathematical  Analysis  of  Visible-Surface  Reconstruction 

l.ct  the  true  distance  from  the  viewer  to  visible  surfaces  be  given  by  the  function  Z(x,y),  where 
x  and  y  are  the  image  coordinates.  I.ow-levcl  visual  processes  generate  a  set  of  noise  corrupted 
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surface  shape  estimates  (i.c..  constraints)  {<',}  which  can  be  expressed  in  the  abstract  notation 

Ci  =  Li[Z(x,  »))  +  £„  (1) 

where  £,  denote  measurement  functionals  of  Z(x.y)  and  (,  denote  associated  measurement  errors. 
Stated  simply,  the  visible-surface  reconstruction  problem  is  to  reconstruct,  as  faithfully  as  possible, 
the  depth  function  Z(x,y)  from  the  available  constraints  {c,}. 

2.1.  The  Ill-Posed  Nature  of  the  Problem  and  Regularization 

The  problem  is  made  nontrivial  by  the  nature  of  the  constraints.  First,  constraints  are  contribt  ted  not 
by  one,  but  by  multiple  specialized  early  visual  processes.  Hence,  slightly  inconsistent  measurements 
provide  1  by  different  prixcsscs  that  happen  to  coincide  will  locally  overdetermine  surface  shape. 
Second,  constraints  are  not  dense,  but  scattered  sparsely  over  the  visual  field.  Therefore,  while 
they  mty  restrict  surface  shape  locally,  they  do  not  determine  it  uniquely  everywhere;  there  remain 
very  m;  ny  feasible  surfaces  that  arc  consistent  with  the  constraints.  Third,  the  measurements  are 
subject  to  errors  and  noise.  High  spatial  frequency  additive  noise,  regardless  how  small  its  (RMS) 
amplitude,  can  locally  perturb  the  surface  (orientation)  radically. 

In  view  of  the  above  three  considerations,  we  cannot  conclude  in  general  that  the  solution 
will  exist,  nor  that  it  will  be  unique,  nor  that  it  will  be  stable  with  respect  to  measurement 
errors,  vlathematical  problems  for  which  die  existence,  uniqueness,  or  stability  of  solutions  cannot 
be  guaranteed  a  priori  arc  said  to  be  ill-posed  (Tikhonov  and  Arsenin,  1977].  Visible  surface 
reconstruction  can  thus  be  characterized  as  a  fundamentally  ill-posed  problem. 

Ill-posed  reconstruction  (or  inverse)  problems  are  the  rule  rather  than  the  exception  :n  early 
vision  [Poggio  and  Torre,  1984],  Ill-posed  problems  cannot  be  solved  in  general,  without  imposing 
some  additional  restrictions  on  possible  solutions.  This  is  the  basis  of  a  number  of  systematic 
approaches,  notably  the  regularization  methods  introduced  by  Tikhonov  and  others  (sec  (Tikhonov 
and  Arsenin,  1977]  and  references  therein).  Duda  and  Hart  [1973,  Sec.  7.4]  mention  a  basic  form 
of  regularization  (essentially  spatial  smoothing)  for  combating  the  effects  of  noise  in  images.  A 
more  sophisticated  class  of  regularization  methods  is  discussed  in  the  context  of  low-level  vision  by 
Poggio  and  Torre  [1984], 

Through  regularization,  ill-posed  problems  can  be  solved  by  reformulating  them  as  variational 
principles  that  arc  effectively  computable.  Unlike  the  original  problems,  the  variational  principle 
formulations  are  well-posed:  i.c.,  it  is  possible  to  guarantee  the  existence,  uniqueness,  and  stability 
of  their  solutions  under  nonrcstrictive  conditions.  Reformulation  proceeds  with  the  introduction  of 
suitable  stabilizing  functionals,  notably  the  class  of  stabilizers  proposed  by  Tikhonov  and  Arsenin 
[1977,  pp.  69-70].  These  stabilizers  can  be  interpreted  as  spline  functionals  that  impose  smoothness 
assumptions  on  die  admissible  solutions  (by  restricting  them  to  Sobolev  spaces  of  smooth  functions).1 
Pragmatically  then,  this  type  of  regularization  is  essentially  equivalent  to  optimal  approximation  by 
generalized  splines  (Tcrzopoulos,  1985a],  We  pursue  the  generalized  spline  approximation  point  of 
view,  since  splines  arc  familiar  and  since  they  suggest  helpful  physical  interpretations. 

2.2.  A  Variational  Principle 

'Ihc  abstract  theory  of  optimal  spline  approximation  is  wcH-devcIopcd  and  a  close  connection  has 
been  established  with  variational  principles  involving  the  constrained  minimization  of  (semi-)  norms 
in  (semi  )  Hilbert  function  spaces  [l  ament,  1972],  Let  f  he  a  linear  space  of  smooth  functions  and 
let  S(n)  be  a  functional  defined  on  V  which  measures  the  (lack  of)  smoothness  of  a  function  in 

1  Generic  smoothness  assumptions  are  generally  the  weakest  (least  committal)  assumptions  that  one  can  make 
about  feasible  solutions  and  still  obtain  well-posed  formuhilions. 
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U.  Furthermore,  let  P  be  a  functional  on  U  which  provides  a  measure  of  the  discrepancy  between 
the  function  and  the  given  constraints.  Consider  the  following  variational  principle: 

VP:  Find  u  €  H  such  that 

£(u)  =  inf  £(v),  (2) 

tifcV 

where  the  energy  functional 

£(v)  =  S(u)  +  P[y).  (3) 

This  variational  principle  will  serve  as  a  formal  statement  of  the  visible-surface  reconstruction 
problem:  The  best  reconstruction  of  the  depth  function  Z(x,y)  from  the  available  constraints  will 
be  given  by  the  solution  u(x,y),  the  smoothest  function  in  the  admissible  space  "H  which  is  most 
compatible  with  the  constraints. 

Before  proceeding  to  specify  the  smoothness  functional  S(v)  and  the  penalty  Junctional  P(v), 
it  should  be  noted  that,  if  the  solution  exists,  it  satisfies  the  necessary  condition  for  the  minimum 
given  by  the  vanishing  of  the  first  variation  6, 

6£(u)  =  6S(u)  +  6P(u)  =  0,  (4), 

which  expresses  the  so  called  Euler- Lagrange  equations. 


2.3.  Generalized  Spline  Functionals 


For  an  appropriate  smoothness  functional  S(v),  we  turn  to  the  multidimensional  splines  studied 
by  Duchon  [1977]  and  Mcinguct  [1979],  generalizations  of  the  classical  univariate  splines  [Ahlbcrg, 
et  al.,  1967].  The  subclass  of  (213)  surface  splines  relevant  to  our  problem  can  be  characterized  as 
members  of  a  suitable  space  of  admissible  functions  v(x,y)  which  minimize  the  functional 


-//,g(7)(=5N  *» 


(5) 


The  positive  integer  m  dictates  the  order  of  the  partial  derivatives  that  occur  in  the  functional, 
which  in  turn  determines  the  order  of  continuity  possessed  by  the  admissible  functions.  The  Bulcr- 
Lagrangc  equation  satisfied  by  the  minimizing  function  u(x,y)  is  an  iterated  version  of  Laplace’s 
equation:  (-l)mAmu  =  0,  where  Au  =  uxx  +  uyy  is  the  Laplacian  of  u. 


Low  order  surface  splines  have  interesting  physical  interpretations  involving  equilibria  of  clastic 
bodies.  Two  special  cases  arc  of  interest.  For  m  =  1  the  functional  reduces  to 

M*  =  //  (*£  +  »})  dxdy,  (6) 

which  is  proportional  to  the  small  deflection  energy  of  a  membrane  (c.g.,  rubber  sheet),  while  for 

m  =  2, 

Ma  =  /  /  ivlx  +  2vly  +  vlv)  dxdV>  (7) 

is  proportional  to  the  small  deflection  bending  energy  of  a  thin  plate  (with  zero  Poisson  ratio) 
[Courant  and  Hilbert,  1953].  Duchon  [1977]  refers  to  the  minimizers  of  [w[|  as  thin  plate  splines. 
Since  thin  plate  splines  arc  the  natural  2D  analogs  of  cubic  splines,  finds  frequent  usage  in 
surface  interpolation  problems  [Schumakcr,  1976].  In  particular,  it  has  been  employed  for  visual 
surface  interpolation  [Grimson,  1983;  Terzopoulos,  1983a]. 


The  physical  interpretations  make  it  clear  that  membrane  splines  offer  a  lower  order  of 
continuity  than  thin  plate  splines.  Since  the  physical  forces  in  the  membrane  arc  due  primarily 
to  its  surface  tension,  it  generates  minimal  area  surfaces.  Although  minimal  area  surfaces  arc 
continuous,  they  need  not  have  continuous  first  partial  derivatives;  i.c.,  they  arc  C’°  surfaces.  For 
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instance  a  sharp  corner  would  result  readily  if  an  ideali/cd  physical  membrane  were  subjected  to 
the  deflecting  force  of  a  knife  edge.  In  contrast  the  restoring  forces  in  a  physical  thin  p  ate  arc 
due  primarily  to  its  flexural  rigidity.  A  thin  plate  would  not  crease  when  deflected  by  a  knile  edge. 
Thin  plate  splines  therefore  maintain  continuity  as  well  as  continuous  first  partial  derivatives;  i.e., 
they  generate  Cl  surfaces. 


2.4.  Controlled  Continuity  and  the  Thin  Plate  Surface  Under  Tension 


Generic  smoothness  assumptions  are  justified  in  pursuing  a  regularization  approach  to  the  visible- 
surface  reconstruction  problem,  inasmuch  as  the  coherence  of  matter  tends  to  give  rise  to  smoothly 
varying  surfaces  relative  to  the  viewing  distance,  over  some  range  of  scales;  however,  smoothness 
assumplions  clearly  do  not  hold  arbitrarily  across  surface  discontinuities,  some  of  which  persist 
across  ell  scales.  1'his  introduces  significant  complications  for  classical  spline  approximation  or 
regulari  ration  methods;  the  continuity  of  spline  functionals  (or  stabilizers)  must  be  controlled  at 
discontiiuitics  in  order  to  preserve  them. 


A  stabilizer  providing  the  necessary  local  continuity  control  can  be  realized  as  a  weighted 
combination  of  generalized  spline  functionals  of  more  than  one  order  m  [Terzopoulos,  1985a].  We 
propose  the  following  smoothness  functional: 


Spriv)  =  ^  J  j  p(*,y)|r(x,y)(v*I  +  2v\y  +  v*y)  +  [1  —  r(x,  y)](v*  +  t^)  j  dxdy,  (8) 

where  f!  denotes  the  image  domain,  and  p(x,y)  and  r(x,y)  are  rcal-.alucd  weighting  Unctions 
whose  range  is  [0, 1],  '1'his  controllcd-continuity  stabilizer  is  a  weighted  convex  combination  of  the 
thin  plate  spline  functional  |tz||  and  membrane  spline  functional  |u|J  integrands.  The  associated 
F.ulcr-Lagrangc  equation  is 


d 

Ox* 


(fXUxx  ) 


+ 


d 

dxdy 


(2puxy)  f  ^  {Puyy) 


d_ 

dx 


(*«.)  -  Qy  (nuy)  =  0, 


(9) 


where  p(x,y)  -  p(x,y)r[x,y)  and  tj(x,y)  -  p{x,y)[l-T(x,y)],  with  natural  (i.e.,  free)  boundary 
conditions.  The  functional  SllT(v)  can  be  thought  of  as  a  thin  plate  surface  under  tension ,  where 
p(x,y)  is  a  spatially  varying  "rigidity”  and  [1  -  t(x,  iy)]  is  the  spatially  varying  “surface  tension.” 
It  generalizes  the  unidimensional  splines  under  tension  of  Schwcikcrt  (see  [Ahlbcrg,  el  ah,  1967]). 


The  local  continuity  properties  of  the  thin  plate  surface  under  tension  functional  can  be 
controlled  at  any  point  (x,  y)  G  U  by  specifying  the  values  of  the  continuity  control  functions 
p(x.y)  and  r[x.y)  at  that  point.  As  r  approaches  1  the  functional  tends  to  a  thin  plate  spline  (a 
(V1  surface)  whereas  towards  the  other  extreme,  0,  the  functional  tends  to  a  membrane  spline  (a 
C’°  surface)  with  intermediate  values  characterizing  a  hybrid  C1  surface  that  blends  the  properties 
of  both  constituent  splines,  p  determines  the  overall  potency  of  the  smoothness  functional. 

Reconstructed  surfaces  must  be  able  to  faithfully  preserve  known  depth  and  orientation 
discontinuities,  while  not  introducing  spurious  discontinuities  at  other  locations.  This  can  be 
accomplished  if  (i)  away  from  known  depth  and  orientation  discontinuities,  the  reconstructed 
surface  possesses  (at  least)  the  Cl  smoothness  of  a  thin  plate,  maintaining  both  continuity  and 
continuous  derivatives,  (ii)  at  known  orientation  discontinuities,  it  exhibits  just  the  C ”  smoothness  of 
a  membrane,  maintaining  continuity  only,  and  (iii)  at  known  depth  discontinuities,  the  smoothness 
functional  is  deactivated  so  that  the  reconstructed  surface  is  free  to  “fracture”  locally.  Hence, 
Spr(v)  will  be  manipulated  as  follows:  At  all  non-discontinuity  points  (x,  y),  p(x,y )  and  r(x,iy) 
should  he  nonzero.  At  orientation  discontinuity  points,  r(x,?/)  is  set  to  zero.  At  depth  discontinuity 
points.  p[j\  y)  is  set  to  zero.  Mechanisms  for  automatically  detecting  discontinuities  hy  computing 
continuity  control  functions  optimally  according  to  local  criteria  arc  considered  in  a  subsequent 
section. 
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2.5.  IV/ialty  Functionals 


Assuming  independently  distributed  measurement  errors  <t  with  zero  means  and  variances  <r2.  the 
optimal  measure  of  incompatibility  is  a  weighted  Kuclidcan  norm  of  the  discrepancy  between  the 
admissible  function  and  the  data  cf. 

W=  5  (i") 

t 

where  the  «,  arc  nonnegative  real-valued  weights  (ideally  a,  is  inversely  proportional  to  o2; 
i.c.,  a,  =  1/A  of)  [Kimcldorf  and  Wahba,  1970],  This  penalty  functional  can  also  be  employed 
(suboptimally)  when  the  above  assumptions  do  not  hold  strictly. 


Appropriate  measurement  functionals  Li  for  surface  reconstruction  may  be  synthesized  from 
generalized  -order  derivatives: 


AM  = 


dkv 


dxJdyk 


'(*• -w.) 


j  =  0,1,...,  k. 


(11) 


k  =  0  yields  simple  evaluation  functionals  A'M(i,y)]  =  tz ( x, ,  2/t ) ,  which  will  be  employed  to 
model  the  local  depth  constraints 

c,  =  v(xi,yi)  +  c,  =  (12) 

The  components  of  the  local  surface  normal  n(:e,,  yt)  =  [vx(x,-,  jtj),t;y(a;,-,j/,), -1],  which  deter¬ 
mine  local  surface  orientation,  can  be  handled  by  the  first  order  ( k  =  1)  derivative  functionals 
£t[t)(x,T/)j  =  vx(xi,y,)  and  Li[v(x,y)]  =  vv(xl,tyl)  and  yield  analogous  expressions  for  the  local 
orientation  constraints: 

c*  =vx(xi,yi)  +  Cj  =  P(ii>y,) 

ci  —vy(Xii  Vi)  +  Cj  —  9(xiiyi)- 

Other  potentially  relevant  functionals  such  as  directional  derivatives  can  be  accommodated  straight¬ 
forwardly  with  the  above  notation. 


It  is  convenient  to  separate  the  various  constraints  into  three  sets;  the  set  i  G  D  of  image 
points  at  which  depth  constraints  rf( Iity.)  occur,  and  the  sets  i  G  P  and  i  €  Q  at  which  orientation 
constraints  P(Xi,y,)  and  q(x,,y,)  occur  respectively.  The  penalty  functional  can  then  be  expressed  as 
a  sum  of  three  components 

p(v)  ZZ\Y adiH^yi)-d(Xi,y.) Y aPi\v*{xi’yi)-P(*<#<)?+\Y a«Wxi>yi)-n*t*< 

ieD  ieP  ieQ 

(14) 

where  the  a,  parameters  are  now  distinguished  as  a^,  aPt,  and  aqi. 


2.6.  A  Physical  Model  for  Visible-Surface  Reconstruction 

The  variational  principle  formulation  of  the  surface  reconstruction  problem  has  an  appealing 
physical  interpretation  which  is  illustrated  in  Fig.  2.  The  thin  plate  surface  under  tension  may 
be  visualized  as  an  clastic  surface,  planar  in  its  natural  state,  whose  clastic  bending  energy  Spr(o) 
stabilizes  surface  shape  so  that  it  varies  smoothly  in  between  constraints  (but  not  at  discontinuities). 
Constraints  deflect  the  surface  according  the  penalty  functional  P(v).  which  can  be  interpreted  as 
the  total  stretching  energy  of  a  set  of  ideal  springs  attached  to  the  constraints.  The  left  part  of 
the  figure  shows  the  clastic  surface  whose  deflection  u(x,y )  at  equilibrium  is  determined  by  an 
infrastructure  of  scattered  depth  constraints.  Ihc  local  depth  estimate  is  encoded  as  the  vortical 
height  of  the  constraint  and  the  tightness  of  each  constraint  is  controlled  by  associated  spring 
stiffness  n.ii.  The  right  part  of  the  figure  illustrates  an  orientation  constraint  coercing  the  local 
surface  normal.  The  spring  stiffness  is  determined  by  the  constraint  parameters  and  rv?1. 
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Figure  2.  The  physical  model.  Thin  plate  surface  under  tension  and  depth  constraints  (left).  Local  influence 
of  an  orientation  constraint  (right). 


2.7.  F xistcnce,  Uniqueness,  and  Stability  of  the  Solution 

Kxistcnce,  uniqueness,  and  stability  of  the  solution  u(x.y)  to  the  variational  principle  VP  are 
guaranteed  when  Spr(v)  =  $pT{v)  +  P(v)  is  a  norm  in  the  admissible  space  X.  Unfortunately, 
generalized  spline  functionals  |w|m  arc  a  priori  only  semi-norms  (of  a  particular  class  of  Sobolev 
spaces).  The  null  spaces  M  of  functions  that  map  to  zero  under  the  semi-norms  are  simply  the 
(M  =  ("‘o' l)  dimensional)  spaces  of  all  polynomials  over  !ft2  of  degree  less  titan  or  equal  to  m  -- 1 
[Meinguct.  1979],  The  penalty  functional  P(v)  can  force  SpT(v)  to  be  a  norm,  however,  if  it  at 
least  constrains  k  to  a  unique  polynomial.  A  possible  set  of  conditions  for  this  to  occur  is  that  the 
£,  include  evaluation  functionals  at  an  .A/ -unisol vent  set  of  points  (i.c.,  a  set  of  M  points  which 
define  a  unique  polynomial  in  the  null  space  of  the  smoothness  functional).  In  particular,  since  the 
maximum  order  of  generalized  splines  in  the  stabilizer  SpT{v)  is  m  =  2,  its  null  space  is  the  space 
of  linear  polynomials.  Thus,  the  following  proposition  can  be  proven  [Terzopoulos,  1984]: 

Proposition.  The  solution  u(x,  y)  will  exist,  be  unique,  and  stable  given  any  one  of  the  following 
minimal  conditions 

(i)  three  noncul inear  depth  constraints, 

(ii)  two  depth  constraints  as  well  as  a  single  p  or  q  constraint, 

(Hi)  a  single  depth  constraint  as  well  as  a  single  p  and  a  single  q  constraint, 

(iv)  a  single  p  and  a  single  q  constraint  with  the  "center  of  gravity"  of  the  surface  fixed. 

These  minimal  conditions  will  hold  in  practice,  due  to  the  large  number  of  constraints  typically 
available  from  early  shape  estimation  processes  (the  fixed  center  of  gravity  condition  can  be  imposed 
when  necessary).  Consequently,  the  visible-surface  reconstruction  problem  may  be  considered  well- 
posed,  hence  effectively  computable  in  general. 

Satisfying  the  conditions  for  a  well-posed  problem  essentially  guarantees  that  a  unique  state  of 
stable  equilibrium  will  exist  for  the  plate/spring  system  (the  minimal  energy  state  £,„■(«)).  In  this 
context,  the  controlled  continuity  assumption  about  surfaces,  as  embodied  by  the  thin  plate  surface 
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under  tension  mode!,  is  physically  nonrcstrictivc  but  nonetheless  powerful  enough  to  guarantee  the 
existence  of  unique  solutions  to  the  variational  principle. 


3.  Discretization 


It  is  extremely  difficult,  if  not  impossible,  to  obtain  an  analytic  solution  to  the  variational  principle 
due  to  the  irregular  occurrence  and  geometry  of  constraints  and  discontinuities.  For  our  purposes, 
the  only  viable  approach  is  to  convert  the  continuous  surface  reconstruction  problem  to  an  equivalent 
discrete  problem  whose  solution  can  be  computed  numerically.  To  this  end,  finite  elements  make 
ideal  local  surface  shape  primitives  for  use  in  visible-surface  representations  [Terzopoulos,  1982, 
1983a].  The  finite  clement  method  [Strang  and  Fix,  1973]  is  a  general,  powerful,  and  mathematically 
rigorous  approximation  technique  which  guides  the  selection  of  appropriate  elements  and  governs 
their  interactions  according  to  the  nature  of  the  variational  principle. 


The  finite  clement  method  offers  substantial  flexibility  in  discretizing  domains  with  irregular 
shaped  boundaries.  Although  the  use  of  irregularly  shaped  elements  to  discretize  such  domains 
may  not  present  a  feasibility  problem  with  regard  to  distributed  biological  mechanisms,  it  makes 
nontrivial  the  mapping  of  elemental  computations  onto  regularly  interconnected  processing  networks 
typically  provided  by  VLSI  technology.  In  this  paper  we  restrict  ourselves  to  regular  finite  elements 
in  order  to  facilitate  such  mappings.  Since  the  goal  is  to  obtain  a  particularly  fine  discretization,  at 
the  resolution  of  the  image,  the  restriction  to  fine  regular  elements  will  not  jeopardize  our  ability 
to  accommodate  the  irregular  occurrence  of  constraints  or  discontinuities. 


3.1.  The  Discrete  Equations 


The  domain  H  is  tessellated  into  square  clement  subdomains  with  sides  of  length  h.  Nodes  are 
located  at  clement  corners  and  shared  by  adjacent  elements.  This  results  in  a  planar  and  uniform 
square  grid  of  nodes  that  is  ideally  suited  to  VLSI  implementation.  The  nodes  arc  naturally 
indexed  by  (i,j)  for  i  =  1  ,...,NX  and  j  =  1  where  Nx  and  Nv  arc  the  number  of 

nodes  along  the  x  and  y  axis  respectively  of  the  (rectangular)  domain  f).  Hie  total  number  of  nodes 
is  N  -  Nx  x  Ny.  The  reconstructed  surface  is  represented  by  an  assembly  of  (nonconforming) 
finite  elements,  each  of  which  is  a  six-point  (full)  quadratic  interpolant  defined  locally  within  its 
particular  subdomain.  The  unknown  displacement  (surface  depth),  at  node  («,  j)  is  denoted  by  die 
variable  •  =  v(ih,jh).  Taken  together,  the  displacement  variables  arc  denoted  by  the  vector 

yh  e  Once  this  vector  is  determined  by  solving  a  discrete  version  of  the  variational  principle, 
the  local  intcrpolants  arc  known  exactly  and,  consequently,  they  explicitly  represent  depth  and 
orientation  everywhere  over  the  surface. 


'Ihc  proposed  square,  quadratic  clement  leads  to  the  following  0 (h2)  formulas  for  the  required 
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partial  derivatives  at  an  arbitrary  node  (i.j)  (Ter/opoulos.  1983a]: 

vu  =*  (v«  »  J.J  "  'V1.;  +  vi  i,j)  ’ 

vyy  =*  (v^  ,i  ~2v^  + 

^  Jfi  (v«'+i,i-»i  v«j+i  ~  v«+i,y  +  V»,y)  !  U**) 

=*•  i  -  <=)  • 

«!  =»  j[  («*,..  -  «?j)  i 

Note  that  the  formulas  arc  finite  difference  expressions.  Their  appearance  is  due  to  the  uniformity 
and  low  order  of  the  clement  and  their  relative  simplicity  will  facilitate  die  calculations  substantially. 

Substituting  the  above  expressions  with  the  constant  approximations  p(x,y)  =>  and 

r(x,y)  =>  t£-  into  (8),  and  noting  that  the  area  of  each  clement  is  hz,  we  obtain  the  discrete 
functional 

S^h)=  v^  +  vtlf/)2 

if  ^ 

+ 2  (v.-n.i+t  -  v?>>+1  -  + <3y 


+ 


+  (v»\y+i  -  2vi,j  +  v«\y-i) 

[‘  -  rij  [(<,.>  -  +  (*i*.  -  v*,)!]  }■ 


(16) 


Although  by  no  means  a  necessity,  it  is  both  natural  (in  view  of  common  image  discretization) 
and  convenient  to  assume  that  die  constraints  coincide  with  nodes  (i,j)  of  the  grid.  Hence,  to 
obtain  a  discrete  expression  for  P(v),  we  collect  the  nodes  at  which  the  various  constraints  occur 
into  three  sets;  die  set  (i.j)  £  D  at  which  depth  constraints  d*;  occur,  and  the  sets  («,  j)  £  P 

and  (i,j)  £  Q  at  which  orientation  constraints  and  q-*y  occur.  Using  symmetric  dilfcrcnce 
approximations  for  the  partial  derivatives  in  (13),  the  discrete  penalty  functional  may  be  written  in 
terms  of  the  nodal  variables  as 


Ihc  energy-minimizing  vector  of  nodal  displacements  uh  satisfies  the  equilibrium  condition 

V£phr(uh)  -  VS^(uh)  +  VP  V)  =  0.  I  >8) 

where  V  is  the  gradient  operator.  Since  the  discrete  functional  ^(u^)  is  a  quadratic  form  in  the 

uthr  the  above  equation  defines  a  linear  system  of  simultaneous  equations  that  arc  satisfied  by  uK 
The  discrete  problem  amounts  to  solving  these  nodal  equations. 
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3.2.  Computational  Molecules 


To  progress  towards  explicit  expressions  for  the  nodal  equations,  we  first  determine  the  partial 
derivatives  of  ^(u*)  and  Ph{ uA)  with  respect  to  an  arbitrary  nodal  variable  u*y.  Letting 


=  piA}/h2  and  =  tfJ1  -  <3)’ 


we  obtain 


_  f  /  *  _  h  h  \  k 

—  |  \ut,y  2ut-i,y  +  Vt-ij 

+  (-2u?+i  j  +  Ks  ~  2u*-u)  & 


dll’?  ■ 


2^1.,  +  at 


,) 


Pi+lJ 


+  (Ui+2,i  - 
+  (Ki  -  2 U.--1J  -  2u!j- 1  +  SuSLjj.j) 

+  (-2u?+lfy  +  2u?,y  +  2u*  -  2u^_,)  /ijj-i 

+  (-2u,-,i+i  +  2u£_|y+J  +  2u£y  -  2u?_w) 

+  (2u-+i,y+i  -  2u*J+1  -  2uJ*+ij  +  2u£y)  /i& 

+  -  2<;-  l  + 

+  (-2ut+i  +  4u»y  -  2u>y_1) 


X*. 


+  (u«,i+2  2Uj  y4  l  +  ut  JJ  /*»,/+ ij 
+  {  (<i  -  li-U  +  (<i  -  U?+V,i) 

+  (ut,j  -  +  (ui,j  -  u«,j+l)  *?«,>  | ’ 


A 


the  discrete  version  of  (9).  Next,  for  (i,  j )  €  D  n  P  fl  Q, 

dPh(uh)  _  h  ..h  „  h  Ah\ 

dut  i  ~  \  d'P  '’i  aji^}) 


+ 


+ 


'  h  n  h  \ 

p»-  lQ  /..ft  _J*  \  _  W-lj  J,  \ 

4ft2  \‘’J  <2'27  2ft  Pil'i) 

(..h  _  h  \,  aPi  hlj  h  ) 

Ah2  \  *>i  +  2 ft  P,+1’J y 


+ 


+ 


'  h  n  h  \ 

4ft2  \  ’’J  ’’J-2/  2ft 

iii  />  _  ,,h  \  .  « J-<  1  nh  \ 

4ft2  l‘'7  *>->'2j  +  2ft 


(19) 


(20) 


(21) 


The  above  expressions  specify  the  nodal  equations  implicitly,  lurch  constituent  term  in  (round) 
parentheses  can  be  represented  graphically  as  a  basic  computational  molecule.  Computational 
molecules  will  be  interpreted  both  as  spatial  representations  of  the  nonzero  coefficients  in  the  nodal 
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equations,  and  as  local  computations  involving  multiplications  and  additions  of  specific  proximal 
nodal  variables.  I'hc  former  interpretation  will  facilitate  the  construction  of  individual  nodal 
equations  given  some  local  structure  of  constraints  and  discontinuities,  while  the  latter  will  lead 
directly  to  local  iterative  algorithms  for  solving  the  resulting  simultaneous  linear  system. 

3.2.1.  Basic  Molecules 


A 


Figure  3.  Plate  molecules  (a)  and  membrane  molecules  (b). 
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3.2.2.  Molecular  Summation  within  Smooth  Regions 

The  formation  of  nodal  equations  within  continuous  regions  can  be  visualized  as  a  prtcess  of 
molecular  summation.  During  molecular  summation,  the  basic  molecules  combine  at  the  central 
node,  coincident  atoms  summing  together. 

When  (i,j)  is  an  interior  node,  away  from  constraints  and  discontinuities,  p^  ■  =  r*.  =  1,  and 
only  the  plate  component  of  (20)  contributes  to  the  expression  for  the  partial  derivative.  Hence, 
the  equilibrium  condition  (18)  reduces  to  the  nodal  equation 

®  ~  ^2  UM  “  ^2  +  U*+l,i  +  U«',y- 1  +  u«,j  +  l) 

+  ^2  +  u?+  i.y-1  +  +  u?+i,y+i)  (22) 

+  h?  +  U‘>2.J  +  U*T/~2  +  UAi+2)  • 

This  equation  can  be  represented  by  the  composite  nodal  molecule  illustrated  in  Fig.  5(a),  which 
results  from  the  summation  of  the  plate  molecules  in  Fig.  3. 


A 


B 


Figure  5.  Interior  node  molecules,  (a)  Away  from  discontinuities,  (b)  At  interior  orientation  discontinuities. 

Note  that  the  computational  molecule  for  the  center  of  the  region  is  a  factor  of  h 2  (due  to  the 
elemental  area)  times  an  order  0(h2)  finite  difference  approximation  for  the  biharmonic  operator 
(Abramowitz.  and  Stcgun,  1%5,  p.  885|,  the  Fulcr-Tag range  equation  associated  with  the  thin  plate 
spline.  This  is  an  expected  consequence  of  the  particular  clement  employed  whicli  yielded  finite 
difference  approximations  for  the  second  partial  derivatives  of  vh. 

If  node  (i,j)  is  a  depth  constraint,  the  first  term  in  (21)  takes,  part  in  the  nodal  equation.  The 
effect  can  be  represented  as  a  summation  of  the  depth  constraint  molecule  and  associated  constraint 
factor  with  the  nodal  molecule  for  (i,j)  shown  in  Fig.  5(a). 

Similarly,  if  (t  -  l,f)  or  (t  +  1.  j)  are  j>  constraints,  or  (i,  j  -  1)  or  (t ,j  +  1)  arc  7  constraints, 
the  other  terms  in  (21)  participate  in  the  nodal  equation.  Again  this  can  be  represented  as  the 
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summation  of  computational  molecules.  Specifically,  the  upper  left  molecule  in  Fig.  4(b)  sums  with 
the  nodal  molecule  if  (*  -  l,j)  G  P.  the  upper  right  only  if  (»  +  1  ,j)  G  P ,  the  lower  left  only  if 
(t ,j  -  1)  G  Q,  and  the  lower  right  only  if  (i,j  +  1)  G  Q. 


3.2.3.  Molecular  Inhibition  at  Discontinuities 


If  (t.  j)  is  a  discontinuity,  cither  or  rj^  or  both  may  be  zero,  thus  nullifying  the  summation 
of  specific  molecules.  This  crucial  influence  of  the  discrete  continuity  control  functions  near  known 
discontinuities  will  be  referred  to  as  molecular  inhibition.  It  was  convenient  for  expressing  (20)  and 
(21)  to  discretize  u(x.y),  p[x.y),  and  r( x.y)  over  the  same  set  of  nodes.  Although  orientation 
discontinuities  can  be  situated  at  these  nodes  (since  u£  •  is  defined  at  an  orientation  discontinuity), 

it  is  better  to  position  depth  discontinuities  on  the  links  half  way  between  nodes  (since  u£  is 
undefined  at  a  depth  discontinuity).  Discretizing  r( x,y)  on  links  docs  not  present  a  problem  in 
practice.  As  a  general  rule,  a  discontinuity  may  inhibit  a  molecule  only  if  it  coincides  with  a 
constituent  atom  or  link. 


First  consider  orientation  discontinuities.  At  an  orientation  discontinuity,  if  -  =  0  and  only 
the  second  component  of  (20)  contributes  to  the  nodal  equation.  In  effect,  the  plate  molecules  arc 
inhibited  and  replaced  by  the  membrane  molecules  of  Fig.  3.  At  an  interior  orientation  discontinuity 
(t.j),  away  from  depth  discontinuities,  all  four  membrane  molecules  superpose  to  yield  the  nodal 
molecule  shown  in  Fig.  5(b),  which  represents  the  nodal  equation 


4u»\,  -  u.-- i,j  ~  "i+ij  ~  1  ~  um  h  =  °- 


The  equation  will  be  recognized  as  - h 2  times  a  standard  finite  difference  equation  for  the  Laplacian 
[Abramowitz  arid  Stcgun,  1965,  p.  885].  It  too  appears  as  a  consequence  of  the  Euler-Lagrange 
equation  associated  with  the  membrane  spline. 


Since  an  orientation  constrain  cannot  meaningfully  coincide  with  an  orientation  discontinuity, 
orientation  discontinuity  nodes  inhibit  orientation  constraint  molecules.  On  the  other  hand,  depth 
constraint  molecules  arc  not  inhibited  by  orientation  discontinuities  since  it  is  perfectly  reasonable 
to  locally  constrain  a  membrane  spline  in  depth. 


Because  smoothness  constraints  arc  unsuitable  at  a  depth  discontinuity  node  (i,j)  (i.e.,  p,-,y  = 
0),  a  nodal  equation  cannot  involve  nodal  variables  separated  by  or  coinciding  with  a  depth 
discontinuity.  Consequently,  depth  discontinuities  inhibit  all  of  the  basic  computational  molecules. 
Fig.  6(a)  illustrates  examples  (disregarding  constraints)  of  nodal  molecules  for  boundary'  nodes 
(marked  as  double  circles)  which  arc  near  depth  discontinuity  nodes  (marked  by  X’s).  Examples  of 
nodal  molecules  at  boundary  orientation  discontinuities  (double  circles)  next  to  depth  discontinuities 
(X's)  arc  shown  in  Fig.  6(b). 


4.  Detection  and  Localization  of  Surface  Discontinuities 


An  important  feature  of  our  framework  for  computing  visible-surface  representations  is  the  uni¬ 
form  treatment  of  constraints  and  discontinuities,  essentially  as  localized  and  independent  surface 
shape  primitives.  This  facilitates  the  parallel  integration  of  discontinuity  information,  along  with 
shape  constraints,  over  the  various  early  shape  estimation  processes.  It  is  convenient  to  think  of 
discontinuity  information  as  being  collected  into  a  discontinuity  map  which  is  in  registration  with 
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the  reconstructed  surface.  Technically,  the  map  comprises  the  nodal  variables  {pf3}  and  {r*;} 
representing  the  discrete  continuity  control  functions. 

Any  early  visual  process  can  participate  in  initializing  the  discontinuity  map  according  to  its 
own  local  hypotheses  about  the  occurrence  of  discontinuities.  In  general,  this  prior  discontinuity 
information  will  be  partly  incomplete  and  inconsistent,  since  it  derives  from  narrowly  specialized  2D 
image  analysis.  In  evolving  a  globally  consistent  surface,  the  visible-surface  reconstruction  process 
performs  a  crucial  task:  it  brings  the  prior  discontinuity  information  into  consonance  with  the  3D 
shape  constraints  collected  from  all  the  early  p rex: esses.  This  raises  the  problem  of  detecting  and 
localizing  both  depth  and  orientation  discontinuities  during  reconstruction.  The  current  section 
investigates  this  problem  for  the  (impoverished)  case  in  which  no  prior  discontinuity  information 
is  available.  We  first  propose  a  straightforward  scheme  which  exploits  the  regularizing  properties 
of  the  surface  model,  then  a  more  sophisticated  approach  that  extends  the  variational  principle  to 
optimally  estimate  discontinuities  according  to  generic  expectations  about  their  local  structure. 

4.1.  Regularization  Based  Discontinuity  Detection 


From  one  perspective,  surface  discontinuity  detection  shares  much  in  common  with  traditional 
approaches  to  image  intensity  edge  detection.  In  particular,  it  is  possible  to  detect  discontinuities 
by  applying  thresholdcd  local  differencing  operations  to  the  reconstructed  surface  which,  like  the 
image,  is  a  regularly  sampled  function.  Because  they  arc  easily  cormptcd  by  image  noise,  however, 
local  edge  operators  such  as  Laplacians  perform  poorly  [Roscnfcld  and  Kak,  1982]  without  a 
smoothing  prefilter,  say  a  Gaussian  [Marr  and  Hildreth,  1980].  Interestingly,  the  thin  plate  surface 
under  tension  performs  the  necessary  smoothing  on  the  sparse  and  noisy  shape  constraints  (standard 
low-pass  filters  such  as  Gaussians  arc  inapplicable  to  sparse  data).  This  regularizing  effect  permits 
the  reliable  computation  of  numerical  derivatives  for  detecting  discontinuities  [Bakhvalov,  1977, 
Sec.  5.4:  Poggio  and  Torre,  1984;  Terzopoulos,  1985a],  In  addition  to  exploiting  the  regularizing 
effect  of  the  thin  plate  surface  under  tension,  the  discontinuity  detection  scheme  described  next 
is  easily  accommodated  within  the  distributed  computational  structure  of  our  framework,  and  it 
permits  relevant  criteria  such  as  psychophysically  measured  limits  on  stcrcofiision  to  impact  on 
discontinuity  detection. 


Consider  the  random  dot  stereogram  in  Fig.  7  which  depicts  a  set  of  planar  surfaces  stacked  in 
depth.  Fig.  8  shows  a  single  continuous  surface  generated  by  the  surface  reconstruction  algorithm 
from  sparse  stereoscopic  disparities  provided  the  Marr-Poggio-Grimson  (MPG)  stereo  algorithm 
[Grimson,  1985].  Fig.  9  dramatizes  a  portion  of  the  reconstructed  surface  in  cross  section  as  it  passes 
across  a  depth  discontinuity.  The  Cl  surface  overshoots  constraints  near  (he  discontinuity  because 
its  smoothness  conflicts  with  the  sudden  change  in  depth.  The  surface  is  clearly  inappropriate  as  a 
final  solution  near  depth  discontinuities,  but  the  local  incompatibility  can  signal  die  occurrence  of 
these  discontinuities. 


Opposing  bending  moments  arc  imparted  to  the  surface  by  the  constraints  on  either  side  of 
the  discontinuity.  The  surface  inflection  (sec  Fig.  9),  where  the  bending  moment  undergoes  a  sign 
change,  localizes  die  depth  discontinuity.  For  a  thin  plate  spline  u(x,y),  the  bending  moment  per 
unit  length  parallel  to  the  x-z  plane  is  proportional  to  -n.„,  while  its  counterpart  parallel  to  the 
y-z  plane  is  proportional  to  ~uyy  [Szilard,  1974],  The  sum  of  orthogonal  bending  moments  gives 
the  total  moment  M  -  -(uxx  +uyy)  =  -Am.  die  negative  l  aplacian  of  the  deflection  function.  It 
can  be  computed  readily  at  a  node  (i,j)  of  the  discrete  surface  using  the  standard  approximation: 
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E# 


Figure  7.  Synthesized  random  dot  stereogram.  When  fused,  the  stereogram  depicts  four  planar  surfaces 
stacked  one  atop  the  other  in  depth. 


The  zero  crossings  of  M  for  the  reconstructed  surface  in  Fig.  8  arc  shown  on  the  left  in 
Fig.  10  as  black  contours.  Most  of  these  correspond  to  weak  inflections  due  to  slight  ripples  in  the 
^constructed  surface.  A  measure  of  significance  is  therefore  needed  to  detect  true  discontinuities 
while  weeding  out  spurious,  weak  inflections.  The  magnitude  of  the  local  depth  gradient  (surface 
tilt)  is  a  suitable  significance  measure  for  depth  discontinuities.  Hence,  an  inflection  point  will 

ft 


be  considered  significant  if  (7  =  jVu| 


z  +  u*  exceeds  a  limit  td  (it  is  more  efficient  to 


use  the  square  of  this  expression  u *  +  uj,  or  even  |ux|  +  ]«y|).  F.mploying  the  usual  discrete 
approximations,  we  obtain 


-  u?  i,i)2  +  -  <>-l)2 


(25) 


The  right  half  of  Fig.  10  shows  the  significant  inflection  points  where  CXJ  >  td  with  td  =  1. 
Adding  these  significant  points  to  the  discontinuity  map  (by  setting  die  associated  to  zero) 
fractures  the  continuous  surface  to  yield  as  a  solution  the  reconstructed  stack  of  surfaces  shown  in 
Fig.  11. 


The  limit  td  must  be  large  enough  so  that  weak  inflection  points  are  rejected  as  possible 
discontinuities,  while  not  so  large  as  to  miss  many  true  depth  discontinuities.  A  possible  criterion 
for  choosing  td  in  applications  to  stcrcopsis  of  opaque  surfaces  is  suggested  by  Panum's  limiting 
case:  i.c.,  when  a  surface  is  tilted  so  much  from  the  viewer  that  it  begins  to  occlude  itself  from  one 
eye,  causing  stcrcopsis  to  fail.  Human  stcreofusion  limits  have  been  measured  psychophysically. 
Using  pairs  of  points  at  different  orientations,  Burt  and  Jules/  (1980]  measured  a  roughly  isotropic 
disparity  gradient  limit  of  approximately  1  between  fusion  and  diplopia.  Interestingly,  this  is  only 
half  the  Panuin  limit. 


It  is  not  inconsistent  with  these  findings  to  use  the  disparity  gradient  limit  t,i  to  detect 
significant  depth  discontinuities  in  conjunction  with  die  isotropic  bending  moments  to  localize 
these  discontinuities.  The  required  local  support  computations  can  be  performed  in  parallel  at  each 
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Figure  8.  Single  surface  reconstructed  from  the  stereogram  of  Fig.  7. 


grid  point  over  the  surface.  Analogously,  significant  orientation  discontinuities  may  be  detected 
when  die  magnitude  of  the  bending  moment  |MU|  of  the  surface  exceeds  a  limit  ta  (points  of 
high  curvature),  and  they  may  be  localized  at  relative  extrema  of  the  bending  moment  (positions 
of  locally  highest  curvature).  The  sign  of  a  bending  moment  extremum  indicates  die  sense  of  the 
orientation  discontinuity;  negative  signals  a  concave  crease,  and  positive,  a  convex  crease.  Curvature 
peaks  were  also  employed  in  a  scheme  for  detecting  surface  orientation  discontinuities  proposed  by 
I^angridgc  (1984). 

4.2.  Discontinuity  Detection  by  Variational  Continuity  Control 

On  the  one  hand,  experimentation  on  natural  data  with  the  regularized  approach  to  discontinuity 
detection  demonstrates  the  feasibility  of  discovering  many  of  the  more  significant  discontinuities 
during  surface  reconstruction  (results  arc  presented  later).  On  the  other  hand,  certain  inherent 
inadequacies  of  this  simple  scheme  ca  often  lead  to  poor  surface  reconstructions.  The  shortcoming 
arc  due  to  a  basic  conflict  caused  by  smoothing.  While  regularization  eliminates  noise,  making 
reasonable  estimation  of  surface  derivatives  possible  in  continuous  regions,  it  tends  to  obscure 
discontinuities  ( I  cr/opoulos,  1985a],  It  can  result  in  poor  detectability  and  localization  of  die  more 
subtle  discontinuities,  a  common  problem  with  smoothing  edge  operators  in  general  [I.cclcrc  and 
Zuckcr,  1984], 
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/*.  significant  inflection 

V  insignificant  inflections 


Figure  9.  Cross-section  of  a  reconstructed  surface  across  a  depth  discontinuity.  The  significant  and  insignificant 
inflectioi  s  of  the  surface  are  indicated. 


The  problem  can  be  resolved  by  exploiting  more  fully  the  controlled  continuity  model  to 
preserve  surface  discontinuities  and,  moreover,  to  incorporate  a  priori  expectations  about  dis¬ 
continuity  structure  into  the  variational  principle  for  surface  reconstruction.  So  augmented,  the 
variational  principle  establishes  a  beneficial  cooperation  between  llic  interpolation  process  which 
smoothly  propagates  shape  information  across  regions,  and  the  complementary  discontinuity  pro¬ 
cess,  which  delimits  these  regions.  Thus  it  optimally  reconstructs  the  piecewise  continuous  turfaccs 
and  discontinuities  simultaneously  to  achieve  the  best  possible  surface  shape. 

As  was  mentioned  in  the  previous  section,  the  smoothness  of  the  thin  plate  surface  under 
tension  is  incompatible  with  any  sudden  transitions  imposed  by  the  scattered  shape  constraints. 
This  implies  that  its  potential  energy  of  deformation  is  generally  greater  at  what  ough,  to  be 
interpreted  as  surface  discontinuities.  Any  local  reduction  in  the  continuity  of  the  surface  'educes 
the  incompatibility  and  locally  reduces  potential  energy.  This  can  be  seen  from  (8):  SpT{v) 
considered  as  a  function  of  (t>.  p.  r),  decreases  as  either  p(x.  y)  or  r(x,  y )  arc  made  zero  over  more 
of  12.  I  his  suggests  that  discontinuities  can  be  discovered  in  the  course  of  solving  the  variational 
principle,  by  allowing  the  surface  to  crease  and  fracture  as  needed  to  reduce  the  total  energy 
below  the  minimum  obtainable  with  a  single  smooth  surface.  ITic  insertion  of  discontinuities  must, 
however,  incur  some  energy  increase,  otherwise  p(x.y)  =  0  everywhere  would  trivially  minimize 
the  energy. 
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Figure  10.  Bending  moment  zero  crossings  (left)  and  detected  depth  discontinuities  (right)  for  the  reconstructed 
surface  in  Fig.  8. 


The  variational  continuity  control  approach  to  detecting  discontinuities  involves  augmenting 
die  original  energy  functional  <f,,r(?>)  -  Sp ,(v)  >  P(v)  with  the  discontinuity  functional  P(p,  t) 
(explained  shortly)  to  obtain  die  new  variational  principle 

Find  u,  p.  and  r  such  that 

£{u,p,t)  =  inf  £(v,p,r),  (2G) 

V,p,T 

where  the  energy  functional 

P(v.p.r)  =  S(o,/?,r)  T  P{v)  I-  P{p,r).  (27) 


The  solutions  u(x p(x.,t/),  and  r(x,y),  satisfy  the  dircc  coupled  F.uler-l.agrange  equations, 
which  express  die  vanishing  of  the  first  variation  with  respect  to  each  independent  function 


fiuS(n.p.  f)  -  0  (llUXX)  +  ~  (2  TiUsy)  +  ~  (llUyy)  ~  ^  fall,)  ~  ~  (’]Uy)  J 


OxOy 

Sl,£(u,p.r)  ■=  0  --T( v2xr  +  2v]y  +  v2yy)  +  (1  7](t>*  +  v2y)  +  £„P(p,r); 

SrS(v.p.T)  -  0  =P  (w*x  +  2v2y  +  v2y)  -  ( V l  +  v2y)  +  6xD(p,t). 

Note  that  the  first  equation  is  identical  to  (9). 


(28) 


Ihc  functional  P(p.r)  maps  the  depdi  and  orientation  discontinuity  configurations  p(: c,y) 
and  r(.r.  y)  into  positive  energies  ( this  is  analogous  to  the  role  of  S,,r(o)  with  respect  to  u).  In  its 
simplest  form,  the  functional  can  increase  monotonically  with  the  total  number  of  discontinuities; 
c.g..  p(p.  t)  -  f  fu /t,/(l  p{s.  </)j  F  ft„[l  -  r(x.y) ]  tlxdy.  where  (itl  and  are  positive  energy 
scaling  parameters  for  the  depth  and  orientation  discontinuity  contributions,  respectively. 


More  interestingly,  significant  advantages  accrue  in  the  detection  of  weak  or  subtle  disconti¬ 
nuities  if  die  functional  can  be  designed  so  as  to  bias  the  solution  according  to  generic  constraints 
about  the  local  structure  of  discontinuities.  Useful  constraints  can,  for  example,  be  based  on 
Gestalt  principles  of  good  continuation  —  discontinuities  tend  to  be  arranged  along  contours,  these 
contours  tend  to  be  continuous,  etc.  This  may  be  accomplished  readily  by  assigning  potential 
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Figure  11.  Reconstructed  surfaces  and  discontinuities 


energies  to  various  local  discontinuity  configurations  on  the  (i.  j)  grid  of  nodes  for  the  discrete 
problem.  I-ncotlings  of  local  edge  configuialions  that  favor  good  continuation  have  been  employed, 
for  instance,  in  relaxation  labeling  curve  enhancement  processes  (/.ticker  ci  <//.,  I977J  and  in  Markov 
random  field  image  restomtion  models  (Genian  and  Geman.  1985], 

In  our  current  implementation,  the  discrete  discontinuity  functional  is  a  weighted  nodal  sum 
of  potential  energy  quanta  !)['  and  ()th}  over  depth  and  orientation  discontinuity  configurations 
respectively: 

;W')  -  f  fio(>i.,(rh)  1.  (29) 

We  employ  a  lion rist it  encoding  which  favors  the  formation  of  continuous  and  smoothly  curving 
contours  by  locally  assigning  higher  energies  to  isolated  discontinuities,  terminations,  sharp  bends, 
junctions,  and  regions,  f  ig.  1?  illustrates  some  of  the  configurations,  and  the  (numeric)  energies 
associated  with  ihem  and  their  rotation. illy  symmetric  counterparts.  Just  as  for  computational 
molecules,  the  circles  denote  nodes  (/./),  while  the  X's  denote  discontinuities  (positions  where  ph 
or  rh  ate  0).  I  he  quanta  m  lug.  12(a)  constitute  Ibi  depth  discontinuities,  which  occur  on 
links  be. ween  nodes  as  explained  previously,  (hey  are  equivalent  to  die  configurations  found  in 
(Geman  and.  Genian.  19S5|.  f  ig.  12(b)  depicts  some  of  the  oiientation  discontinuity  configurations 
encoded  b\  ()[  ..  Orientation  discontinuities  coincide  with  nodes. 
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ininimi/dtion.  The  nonconvcxity  of  ihc  energy  landscape  makes  this  a  much  more  difficult  problem 
to  solve  numerically.  In  the  discontinuity  detection  experiments  to  be  presented  in  See.  fi.5,  we 
propose  a  strategy  for  efficiently  obtaining  good,  though  not  necessarily  optimal  solutions. 


5.  Overview  of  the  Multiresolution  Surface  Reconstruction  Algorithm 


Application  of  the  finite  element  yields  the  discrete  problem  of  solving  a  linear  system  of  simul¬ 
taneous  equations.  This  system  has  computationally  desirable  properties;  i.e.,  its  matrix  is  sparse, 
banded,  symmetric,  as  well  as  positive  definite  (for  fixed  p(x,y)  and  tau(x,y))  when  the  available 
constraints  satisfy  the  conditions  for  a  well-posed  problem.  Ihc  sparseness  of  the  matrix,  a  direct 
consequence  of  the  local  support  of  die  finite  element,  is  evident  from  the  nodal  equation  of  an 
interior  node;  rows  associated  with  interior  nodes  have  only  13  nonzero  entries,  while  nodes  at  and 
near  discontinuities  have  even  fewer.  The  N  x  N  matrix  however,  tends  to  be  extremely  large  in 
practice,  since  the  number  of  pixels  N  in  a  typical  image  can  range  from  104  to  106  or  greater. 
This  combination  of  properties  suggests  the  application  of  iterative  techniques  such  as  (parallel) 
Jacobi  or  (sequential)  Gatiss-Scidcl  relaxation  methods  [Hageman  and  Young,  1981].  Relaxation 
methods  lead  to  distributed  algorithms,  and  the  parallel  variants  may  be  implemented  concurrently 
on  networks  of  many  simple,  locally-interconnected  processors. 


5.1.  Nodal  Relaxation  Computations 


A  local-support  nodal  relaxation  computation  can  be  obtained  at  node  (t ,j)  by  expressing  u£y 
in  terms  of  the  remaining  vat  .ablcs  in  the  nodal  equation  determined  by  the  local  structure  of 
constraints  and  discontinuities.  The  nodal  relaxation  computation  may  be  constructed  automatically 
by  applying  our  simple  rules  governing  the  summation  of  basic  computational  molecules: 

(i)  Plate,  depth  constraint,  and  orientation  constraint  molecules  sum  at  interior  (non-discontinuity) 
nodes. 

(ii)  Membrane  and  depth  constraint  molecules  sum  at  orientation  discontinuity  nodes. 

(iii)  Orientation  discontinuities  inhibit  plate  and  orientation  constraint  molecules. 

(iv)  Depth  discontinuities  inhibit  all  basic  molecules. 

For  instance,  at  a  depth  constraint  node  away  from  discontinuities,  the  Gauss-Scidcl  relaxation 
computation  becomes 

8_ 

/»2 

_  1  /\]("  D)  +u(n+1>  4-  n(n)  +u(n)  1 

/,2  Vu*  cj-i  +  u>n,j  -i  +  ut-i,y+i  +  un 


(ml) 

U  •  — 

t,j 


1 


!)  ,„(")  (  ,i(n  liln)  ) 

+  u«ll,y  +  ui,y  1  +  ui,y+lj 


(30) 


1  (iSni  4-  1|M  4-  ill”'1  4-  i/n)  d •• 

2,j  +  U«>2,j  +  ij  2  +  U*  4  2  y  + 


where  we  have  suppressed  the  discretization  superscript  h  and  instead  introduced  the  bracketed 
iteration  indices.  At  an  unconstrained  depth  discontinuity,  we  obtain 


•  t) 


=  i/u(nni 

4  V  1  l-J 


l)  4.  „(») 


-i,y  ^  y”«  i,}  ^  U«,y- i 

Note  that  the  nodal  relaxation  computations  do  not  change  from  one  iteration  to  the  next,  so  long 
as  the  influencing  constraints  or  discontinuities  remain  unperturbed. 


(«+i) 


u(n)  ) 

u*,y  1 1/  * 


(31) 
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S.2.  Multircsolution  Relaxation 

A  serious  problem  with  iterative  techniques,  in  general,  is  their  slow  convergence  rates  for  large 
problems.  This  inherent  inefficiency  is  due  to  the  fact  that  information  must  propagate  incrementally 
across  large  representations  from  nodes  to  their  near  neighbors  in  accordance  with  the  nodal 
relaxation  formulas.2  We  have  developed  highly  efficient  iterative  algorithms  that  overcome  this 
problem  for  surface  reconstruction  [I'erzopoulos,  1982,  1983a]  as  well  as  for  certain  other  visual 
problems  [ I’erzopoulos.  1984],  These  algorithms  achieve  efficiency  by  exploiting  multircsolution 
relaxation  methods  [Fedorenko,  1961;  Brandt,  1977,  Hackbusch  and  Troltenbcrg,  1982], 

Briefly,  the  multiresolution  surface  reconstruction  algorithm  features  (i)  multiple  representations 
of  surface  shape  over  a  range  of  spatial  resolutions,  (ii)  local,  iterative  (relaxation)  processes 
that  propagate  smoothness  constraints  within  each  representational  level,  (iii)  local  coarse- to-fine 
(prolongation)  processes  that  allow  coarser  representations  to  constrain  finer  ones,  (iv)  finc-to-coarse 
(restriction)  processes  that  allow  finer  representations  to  constrain  and  improve  the  accuracy  coarser 
ones,  and  (v)  a  multilevel  coordination  strategy  that  enables  the  hierarchy  of  representations  and 
component  processes  to  cooperate  towards  increasing  the  computational  efficiency,  usually  by  orders 
of  magnitude. 

Fig.  13  depicts  the  structure  of  the  algorithm  schematically.  In  this  particular  case,  only  three 
levels  arc  shown.  Note  the  2: 1  resolution  reduction  between  adjacent  levels.  Not  only  docs  this 
ratio  simplify  the  component  processes  considerably,  but  it  is  also  nearly  optimal  with  regard  to 
total  computation  to  convergence  (this  is  conveniently  measured  in  machine  independent  work 
units ,  where  a  work  unit  is  the  amount  of  computation  required  for  a  relaxation  iteration  on  the 
finest  level)  [Brandt,  1977].  The  diagram  illustrates  the  intralcvcl  relaxation  processes,  as  well  as  the 
fine-to-coarsc  restriction  and  coarsc-to-finc  prolongation  processes  that  communicate  between  levels. 
The  figure  shows  synthetically  generated  scattered  orientation  and  depth  constraints  consistent  with 
a  hemispherical  surface.  The  algorithm  reconstructs  a  dense  representation  of  surface  at  three 
resolutions.  The  sparse  information  at  any  particular  scale  can  be  thought  of  as  a  set  of  constraints 
which  defines  a  discrete  surface  approximation  problem  at  that  level.  It  is  natural  then  to  view  the 
multiresolution  surface  reconstruction  algorithm  as  iteratively  solving  a  coupled  hierarchy  of  discrete 
surface  reconstruction  problems.3  For  a  detailed  description  of  the  algorithm  see  [Terzopoulos, 
1982,  1983a]. 


6.  Experimental  Analysis  of  the  Algorithm 

The  multircsolution  visible-surface  reconstruction  algorithm  was  tested  on  a  variety  of  data  sets 
including  synthetic  data,  structured  light  (laser)  range  data,  automated  stercopsis  and  photometric 
stereo  data  from  natural  images,  and  digital  terrain  model  data.  Some  results  arc  presented 
in  this  section  (for  further  details  and  examples,  sec  [Terzopoulos,  1984]).  In  all  the  examples 

2  It  is  possible  to  accelerate  the  basic  relaxation  methods  so  that  fewer  iterations  are  required.  However,  prac¬ 
tical  accelerated  methods  such  as  the  conjugate  gradient  method,  successive  overrelaxation,  and  Chcbyshev 
semi-iteration  use  global  procedures  to  determine  the  acceleration  parameters.  In  parallel  implementation, 
the  greater  complexity  of  the  globally  accelerated  methods  and,  even  more  importantly,  the  communications 
costs  of  perfonning  the  global  operations  nullifies  any  potential  gains. 

3  A  recursive  multilevel  coordination  strategy  was  employed  in  the  experiments  described  next  The  recursive 
strategy  activates  only  a  single  level  at  any  one  time.  We  have  recently  developed  a  concurrent  strategy 
based  on  a  multilevel  variational  principle  [  Terzopoulos,  1985b].  Concurrent  coordination  maintains  all 
levels  active  simultaneously,  thus  achieving  full  parallelism. 
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presented,  the  intralevcl  process  was  Gauss-Scidcl  relaxation  and  the  algorithm  was  started  from 
zero  initial  approximations  on  all  the  levels.  The  border  nodes  on  each  grid  were  preset  as  depth 
discontinuity  nodes  to  introduce  natural  boundary  conditions  which  free  the  reconstructed  surface 
on  the  boundary  of  ft. 

6.1.  Synthetic  Data 

The  first  two  examples  involve  randomly  placed  depth  constraints.  The  left  half  of  Fig.  14 
shows  15%-dcnsity  constraints  at  three  resolutions,  'lhcsc  constraints  were  obtained  by  sampling 
a  hemisphere  whose  z  values  were  multiplied  by  a  radial  sinusoid.  'Ihe  nodes  outside  the  circular 
region  occupied  by  the  constraints  were  specified  as  depth  discontinuities.  'ITtc  reconstructed  surface 
representation  is  shown  on  the  right  half  of  Fig.  14.  In  Fig.  15,  the  15%-dcnsity  depth  constraints 
shown  on  the  left  arc  samples  of  a  stacked  set  of  planar  surfaces  at  three  resolutions.  In  this 
example,  depth  discontinuities  were  placed  along  the  circular  arcs  bounding  the  planes,  and  along 
the  outer  edges  of  the  grids.  Ihe  rcconstaictcd  surface  representation  is  shown  on  the  right  half  of 
Fig.  15.  This  example  indicates  that  discontinuities  can  be  placed  along  arbitrary  contours  within  ft 
to  prevent  surface  shape  from  being  degraded  by  unwanted  smoothing  over  sharp  depth  changes. 

The  next  examples  involve  reconstructions  from  orientation  constraints.  'Ihe  left  half  of  Fig.  16 
shows  in  perspective  a  set  of  orientation  constraints  over  a  square  region.  On  each  of  three  scales, 
the  region  is  divided  into  four  quadrants  each  containing  constant  orientation  constraints,  and  the 
nodes  along  their  boundaries  arc  preset  as  orientation  discontinuities.  Ihe  surfaces  reconstructed 
by  the  three-level  algorithm  are  shown  on  the  right  half  of  the  figure.  Since  absolute  depth  cannot 
be  determined  solely  from  orientation  constraints,  a  relative  depth  reconstrucdon  results,  with  the 
center  of  gravity  of  the  resulting  pyramidal  surface  resting  near  the  x-y  plane. 

Ihe  left  half  of  Fig.  17  shows  30%-dcnsity  scattered  orientadon  constraints  consistent  with  a 
hemispherical  surface  at  three  resolutions.  The  reconstructed  surface  representation  is  shown  on 
the  right.  All  nodes  outside  the  hemispherical  surface  patch  were  specified  as  depth  discontinuities. 
Again,  the  center  of  gravity  of  the  surface  rests  near  the  x-y  plane. 

The  next  examples  demonstrate  die  integration  of  both  depth  and  orientation  constraints.  The 
left  half  of  Fig.  18  shows  15%-dcnsity  depth  constraints  consistent  with  a  hemispherical  surface  at 
three  resolutions.  On  the  right  arc  15%-dcnsity  orientation  constraints  consistent  with  the  same 
surface.  Nodes  outside  the  surface  have  been  specified  as  depth  discontinuities.  I  he  reconstructed 
surface  is  shown  in  Fig.  19.  Whereas  in  the  previous  example  (Fig.  17)  only  relative  depth  can  be 
determined  for  lack  of  any  depth  constraints,  in  the  present  example  the  additional  depth  constraints 
enable  the  absolute  depth  of  the  surface  to  be  determined  at  all  points,  hence  the  surface  is  “raised” 
to  the  correct  height  above  the  base  plane.  In  addition,  note  that  (10%)  uniformly  distributed 
noise  has  been  added  to  the  constraint  values.  With  the  given  constraint  parameters,  the  surface  is 
slightly  bumpy  on  the  finest  level.  This  can  be  reduced  by  decreasing  the  constraint  parameters,  in 
effect  loosening  the  springs  of  the  physical  model. 

6.2.  Structured  Light  Data 

Ihe  multircsolution  algorithm  was  applied  to  (lie  reconstruction  of  several  objects  from  raw 
range  data  supplied  by  a  laser  rangefinder  constructed  by  P.  Brou  at  MI  T.  I  he  scan  resolution 
in  the  y  direction  is  half  that  in  the  x  direction.  A  four-level  surface  reconstruction  algorithm 
was  employed  in  die  examples.  'Ihe  data  was  introduced  as  depth  constraints  at  the  finest  level 
and  transferred  to  the  coarser  levels  by  successive  2x2  averaging  between  levels.  To  expediently 
segment  the  objects  from  the  background,  values  smaller  dian  a  threshold  were  treated  as  depth 
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Figure  16.  Rcconsiruclion  of  a  pyramidal  surface  with  orientation  disconlimiities  from  oricnuition  constraints. 
(Grid  dimensions:  N’j'  =  N’f  --  17.  Af'11  =  =  33.  W'13  --  -  G5.  Grid  spacings:  /i,  =  0.4, 

/i2  -  0.2.  /i3  -  0.1.  Constraint  parameters:  =  oj'  =  4.0//»;.  Computation:  19.5  work  uniLs.) 
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Figure  20.  Reconsiruciion  of  a  lightbulb  from  range  data.  (Finest  grid  dimensions:  N*'  x  N*'  -  257  x  281. 
Grid  spacings:  J»i  =  0.8.  h2  =  0.4,  h2  -  0.2.  and  64  -  0.1.  Constraint  parameters:  ajh>  =  0.2/hj. 
Computation:  9.78  work  units.) 


discontinuities.  Fig.  20  shows  the  reconstructed  surface  of  a  lightbulb.  The  algorithm  smoothes  the 
noise  in  the  data  and  reconstructs  the  missing  points. 

6.3.  Natural  Image  Data 

In  this  section,  we  apply  the  mullircsolution  surface  reconstruction  algorithm  to  depth  data  orig¬ 
inating  from  natural  images.  'Ihc  examples  involve  photometric  stereo,  and  two  binocular  stereo 
algorithms  applied  to  tcirain  stcrcopairs. 

6,3.1.  Photometric  Stereo  Data 

Photometric  stereo  is  a  technique  that  uses  multiple  (usually  3)  images  of  a  scene  from  the  same 
viewpoint,  but  with  dilfering  illumination  [Woodham,  198 1 j.  Assuming  that  the  surface  material 
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Figure  21.  Image  of  a  maile  white  torus. 


is  known  and  that  the  viewer  and  lightsourccs  arc  far  from  the  object,  the  method  determines 
the  surface  orientation  from  the  image  irradiancc.  Our  surface  reconstruction  algorithm  provides 
a  noise  resistant  technique  for  computing  depth  from  the  surface  orientation  data  provided  by 
photometric  stereo.  We  demonstrate  this  with  fire  image  of  tom s  in  Fig.  21.  The  photometric  stereo 
data  (generated  by  a  system  implemented  at  MIT  by  K.  Ikcuchi)  was  introduced  as  orientation 
constraints  on  a  two-level  algorithm.  Aside  from  sporadic  missing  data,  the  constraints  on  the 
coarse  level  arc  dense,  whereas  only  every  other  node  on  the  fine  level  is  a  constraint.  Fig.  22 
shows  the  orientation  data  and  the  reconstructed  torus. 

Our  method  for  reconstructing  surfaces  from  scattered  orientation  constraints  can  be  compared 
to  a  variational  scheme  for  obtaining  relative  depth  from  dense  surface  gradient  infoimation  reported 
by  Horn  and  Brooks  (1985],  Their  proposed  least  squares  integral  f  f(vx-  p)2  +  (vy  -  q)2  dxdy 
will  be  recognized  as  being  a  continuous  version  of  the  orientation  constraint  penalty  functional.  By 
virtue  of  the  additional  smoothness  functional  however,  our  surface  reconstruction  algorithm 

can  deal  with  orientation  constraints  that  arc  scattered.  It  also  can  integrate  depth  constraints  from 
other  sources  to  arrive  at  absolute  surface  depth. 

6.3.2.  Correlation  Based  Stereo  Data 

At  the  top  of  Fig.  23  is  a  stercopair  on  which  Kass’s  [1983)  correlation  based  stereo  algorithm 
was  run.  The  output  of  the  stereo  algorithm  is  shown  on  the  lower  left,  with  brightness  proportional 
to  disparity .  The  algorithm  has  failed  to  produce  a  match  in  the  neutral  grey  patches,  so  disparity 
is  unknown  in  these  areas.  To  apply  the  multircsolution  algorithm,  the  disparity  data  on  the  finest 
level  were  reduced  by  factors  of  two,  through  averaging,  to  three  coarser  levels.  Relatively  small 
constraint  parameter  values  were  chosen  in  order  to  counteract  die  potentially  detrimental  effects 
of  false  matches  and  noise  in  the  disparity  data.  The  reconstructions  on  the  three  coarsest  levels 
arc  shown  as  3f)  plots  in  Fig.  24  (the  finest  level  was  too  dense  to  represent  this  way).  Fig.  25 
shows  isoclcvation  contour  maps  of  the  solution  on  all  levels. 
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Figure  22.  Reconstruction  of  the  torus  (right)  from  the  orientation  constraints  provided  by  photometric  stereo 
(left).  (Grid  dimensions  JV* 1  =  N'f'  =  51  ;md  N^‘  -  N ^  =  101.  Constraint  parameters:  ah>  =  4.0 /hj. 
Computation:  52.0  work  units.) 


6.3.3.  Feature  Based  Stereo  Data 


The  next  example  involves  disparity  constraints  generated  by  the  MPG  stereo  algorithm 
[Grimson,  1985J.  A  three-channel  version  of  the  stereo  algorithm  was  run  on  the  stcrcopair  at 
the  top  of  Fig.  26.  The  output  of  the  stereo  algorithm  is  shown  on  the  lower  part  of  the  figure. 
Disparity  information  is  provided  only  along  zero  crossing  contours  at  the  three  finest  scales.  In 
the  figure,  the  darkness  along  contours  is  proportional  to  disparity.  This  disparity  data  was  input  to 
a  four-level  surface  reconstruction  algorithm.  The  constraints  on  the  coarsest  level  were  derived  by 
averaging  die  constraints  from  the  next  finer  level.  The  reconstructions  on  the  three  coarsest  levels 
arc  shown  as  31)  plots  in  I  'ig.  27.  Fig.  28  shows  isocicvalion  contour  maps  of  the  solution  on  all 
levels. 


6.4.  Digital  Terrain  Map  Data 
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Figure  23.  Natural  terrain  stercopnir  (top)  and  output  of  Kass'  stereo  algorithm  (bottom).  The  images  were 
25G  x  25G  pixels,  quantized  to  256  levels  (provided  by  the  US  Defense  Mapping  Agency). 


Figure  24.  Reconstruction  of  terrain  in  Fig.  23.  (Grid  dimensions:  N*'  -  N*'  -  33,  N*2  =  N**  ---  G5, 
=  JV*-1  -  129.  N1;*  =  iV'N  -  257.  Grid  spacings:  -  0.8.  h2  =  0.4.  =  0.2,  fc,  ~  0.1. 

Constraint  parameters:  or/1'  O.OJ/iiJ.  Compulation:  29.0  work  units.) 
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Figure  25.  Isoelevation  contour  maps  of  the  reconstructed  terrain  in  Fig.  24. 


A  four-level  surface  reconstruction  algorithm  was  applied  to  contoured  terrain  elevation  data. 
A  contour  map  of  die  Black  River  Gorges  (published  by  die  UK  Ministry  of  Defense)  was 
digitized  manually  on  a  digitizing  tablet  by  J.  Mahoney.  The  25G  x  25G  digital  contour  array  is 
shown  at  die  top  of  l-’ig.  29.  The  constraints  input  to  die  algorithm  are  shown  at  die  bottom 
of  the  figure.  The  elevation  of  die  contours  is  proportional  to  brightness.  Local  averaging  was 
used  to  derive  the  constraints  on  the  coarser  grids  from  diose  on  the  finest  grid.  The  terrain 
reconstructions  on  the  three  coarsest  levels  arc  shown  as  31)  plots  in  Tig.  30  (the  finest  level  is  too 
dense  to  represent  diis  way).  Fig.  31  shows  isoelevation  contour  plots  of  the  reconstructed  terrain 
on  all  levels.  The  reconstructed  contours  on  the  finest  level  can  be  compared  subjectively  with 
die  digitized  contours  in  Fig.  29.  but  note  the  reconstructed  contours  depict  elevations  half  way 
between  the  original  constraint  contours  for  an  unbiased  comparison.  The  reconstructed  contours 
arc  somewhat  smoother  than  die  (predigitized)  contours  in  die  original  map  —  the  jaggedness 
introduced  by  manual  digitization  has  been  reduced.  The  extent  of  t’  nothing  can  be  regulated 
by  adjusting  the  constraint  parameters.  Shaded  image  renditions  of  die  reconstructed  terrain  using 
reflectance  map  techniques  for  hill  shading  [Morn.  198 1 J  arc  shown  at  the  bottom  of  Kig.  31.  Terrain 
reconstructions  using  the  thin  plate  surface  under  tension  model  were  compared  to  reconstructions 
using  die  simpler  membrane  spline  model  (Laplacian  smoothing).  'Hie  former  gives  good  results, 
whereas  the  latter  generally  suffers  from  insufficient  smoothness  and  produces  flat  spots  across 
terrain  peaks  [Terzopoulos,  1984J  (sec  also  [Bolondi  el  al.,  1976]). 
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Figure  26.  Natural  terrain  stereopair  (top)  and  output  of  the  MPG  stereo  algorithm  (bottom).  The  images 
were  512  x  512  pixels,  quantized  to  256  levels  (provided  by  the  US  Army  Engineer  Topographic  Labs). 
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Figure  27.  Reconstruction  of  data  in  Fig.  26  on  the  three  coarsest  levels.  (Grid  dimensions:  TV'*1  =  N'j' 
Nh,  _  Nh,  _  C5  ^t.3  =  Nh,  -  |29,  JV**  -  =  257.  Grid  spacing::  hx  =  0.8,  /i2  =  0.4,  h3 

=  0.1.  Constraint  parameters:  a,;'1'  =  0.01//iJ.  Computation:  31.0  work  units.) 


I 
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Figure  28.  Isoelevation  contour  maps  of  the  reconstructed  terrain  in  Fig.  27. 


6.5.  Discontinuity  Detection  Experiments 

The  foregoing  examples  have  shown  that  the  surface  reconstruction  algorithm  can  handle  disconti¬ 
nuities  that  arc  prcspccificd.  In  the  next  two  sections,  we  present  examples  involving  the  automatic 
detection  of  discontinuities. 

6.5.1.  The  Regularization  Approach 

The  aerial  view  stcrcopair  of  Fig.  32  was  input  to  the  MPG  stereo  algorithm  which  generated 
the  sparse  disparity  map  shown  on  the  top  of  Fig.  33.  The  finest  level  dense  disparity  map  generated 
by  a  four-level  surface  reconstruction  algorithm  is  shown  at  the  iower  left  of  the  figure.  Darkness 
is  proportional  to  disparity.  The  discontinuities  found  from  this  disparity  map,  using  a  disparity 
limit  >  tj  -  1  arc  shown  at  the  lower  right  as  white  contours  After  the  detected  points  arc 
added  to  the  discontinuity  map,  the  surface  reconstruction  algorithm  continues  iterating  from  the 
tentative  approximation  on  the  left.  The  amount  of  additional  computation  required  is  relatively 
small,  since  the  tentative  surface  is  a  fairly  good  approximation  in  most  places.  At  convergence,  the 
reconstructed  surface  has  fractured  along  the  contours  to  give  the  solution  on  the  right.  Portions 
of  the  main  discontinuities  around  the  buildings  have  been  found,  but  contours  arc  broken  and 
shifted. 

The  next  example  involves  the  synthesized  random  dot  stereogram  in  Fig.  7.  The  depth 
constraints  generated  by  a  thrcc-channel  version  of  the  MPG  stereo  algorithm  arc  shown  in  Fig.  34 


Figure  29.  Digitized  contour  data  (top)  and  constraints  (bottom).  The  patch  to  the  lower  right  represents  a 
lake. 
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Figure  30.  Reconstruction  of  data  in  Fig.  29.  The  terrain  reconstruction  on  the  three  coarsest  levels  is 
represented  as  3D  surface  plots.  (Grid  Dimensions:  N*‘  =  N*'  =  33.  N'f1  =  N£J  =  C5,  JV*’  =  N*3 
129.  IV**  =  N*4  -  257.  Grid  spacings:  =  0.8,  h3  --  0.4.  h3  =  0.2,  /»<  =  0.1.  Constraint  parameters: 


Terzopoulos 


Computing  Visible-Surface  Representations 


46 


Figure  31.  Isoclevation  contour  map  (lop)  of  the  data  in  Fig.  30  and  shaded  representations  of  the  reconstructed 
terrain  (bottom). 
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Figure  32.  Aerial  view  of  a  hospital  complex.  The  slereopair  was  provided  by  the  UBC  Faculty  of  Forestry. 
Images  are  320  x  320  pixels. 


(the  finest  level  dimensions  are  320  x  320).  The  constraints  on  the  coarsest  level  were  obtained  by 
averaging  those  on  the  next  finer  level.  Fig.  35  shows  the  smooth  disparity  maps  initially  computed 
by  a  four-level  surface  reconstruction  algorithm.  Fig.  36  shows  the  discontinuities  detected  from 
these  maps  with  tj  =  1.  The  discontinuities  have  been  superimposed  onto  the  final  disparity  maps 
in  Fig.  37.  Better  performance  is  observed  in  this  case  due  to  the  simpler  surface  structure,  but  the 
contours,  while  mostly  intact,  arc  quite  ragged. 

In  general,  not  detecting  true  discontinuities  affects  surface  shape  more  advciscly  over  larger 
regions  titan  introducing  some  spurious  ones  within  a  continuous  surface.  Discontinuity  points 
are  missed  by  the  thresholding  operation,  and  no  adjustment  of  the  global  limit  can  be  expected 
to  produce  perfect  results.  Note,  however,  that  the  surface  reconstruction  algorithm  docs  not 
break  down.  Rather,  the  reconstructed  surface  degrades  as  it  “leaks”  through  the  gaps.  The 
discontinuity  detection  procedure  may  be  improved  by  allowing  the  disparity  limit  to  vary  spatially, 
or  by  modifying  it  during  multiple  passes.  On  the  first  pass,  surface  shape  is  poorest,  so  a  fairly 
conservative  limit  should  be  set  to  reduce  the  number  of  false  detections.  Conservative  limits  fail 
to  detect  many  discontinuities,  but  as  more  discontinuities  are  identified,  surface  shape  improves 
and  limits  can  be  lowered  in  subsequent  passes  to  find  the  less  prominent  discontinuities. 

6.5.2.  The  Variational  Continuity  Control  Approach 


A  multipass  scheme  is  also  employed  in  the  variational  continuity  control  approach  to  efficiently 
obtain  good  solutions.  An  example  will  be  used  to  explain  the  strategy.  Fig.  38  shows  depth 
constraints  randomly  sampled  from  a  set  of  sloping  planes  that  form  discontinuities  along  their 


Figure  33.  Discontinuities  in  the  aerial  stereogram.  Disparity  contours  generated  by  stereo  algorithm  (top), 
full  disparity  map  generated  by  the  surface  reconstruction  algorithm  at  the  finest  level  (lower  left),  and  delected 
discontinuities  superimposed  on  the  disparity  map  (lower  right). 


I 
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Figure  34.  Depth  constraints  for  the  random  dot  stereogram. 


extremities.  A  single  continuous  surface  can  be  rcconstmctcd  from  these  constraints  as  shown,  but 
it  smooths  over  the  depth  discontinuities  and  rounds  out  the  orientation  discontinuities. 

I  hc  algorithm  finds  both  kinds  of  discontinuities  and  reconstructs  a  surface  which  preserves 
them.  When  die  surface  smooths  through  a  depth  discontinuity,  two  spurious  regions  of  high 
curvature  border  the  discontinuity,  'these  spurious  regions  can  easily  be  mistaken  for  orientation 
discontinuities.  To  avoid  this  unwanted  interaction  which  can  substantially  slow  down  the  optimiza¬ 
tion  process,  the  algorithm  postpones  the  orientation  discontinuity  detection  phase  until  all  depth 
discontinuities  have  been  found.  The  surface  evolves  in  several  steps  over  which  the  parameters 
and  ft*  in  (29)  are  modified.  Hach  step  consists  of  first  flipping  the  value  of  the  continuity 
control  parameter  (p*}  or  t*})  from  0  to  1  or  conversely,  if  this  lowers  the  energy  (27),  and  then 
miming  the  reconstruction  algorithm  to  convergence  (which  always  results  in  equilibrium,  since  the 
variational  principle  is  convex  for  fixed  p*j  and  r*-). 

For  depth  discontinuities,  ft*  is  initially  set  to  a  high  value  that  heavily  penalizes  their  insertion, 
then  lowered  in  steps.  This  strategy  of  least  commitment  finds  the  prominent  discontinuities  earliest, 
improving  the  surface  as  it  docs  so,  and  leaves  the  more  subtle  ones  for  later.  It  results  in  the 
(lipping  of  relatively  few  variables  in  each  stage,  hence  the  solution  is  obtained  efficiently,  beginning 
with  the  continuous  surface  Fig.  39(a),  Fig.  39(b-d)  illustrates  the  steps  of  llic  evolving  discontinuity 
detection  process,  during  which  discontinuities  arc  determined  with  increasing  accuracy  as  ft)  is 


lowered.  The  energy  can  be  lowered  further  still  if  /I*  is  then  increased  slightly  to  eliminate 
spurious  discontinuities  in  I  'ig.  39(d).  Note  that  since  die  surfaces  have  now  separated,  a  very  large 
increase  would  be  needed  to  flip  a  true  discontinuity  point  (a  hysteresis  effect).  The  improved 
surface  in  fig.  39(e)  results.  Next,  the  orientation  discontinuity  detection  phase  is  activated  and  it 
runs  in  die  same  way,  but  modifies  In  diis  example,  the  orientation  discontinuities  arc  found 
in  only  one  step. 

fig.  39(f)  shows  (lie  final  solution.  The  depth  and  orientation  discontinuities  have  been  made 
explicit  and  arc  preserved  by  the  reconstructed  surface.  Incidentally,  the  global  optimum  of  the 
variational  principle  has  been  found  in  this  example:  however,  this  procedure  can  generally  be 
expected  to  yield  good,  though  not  necessarily  optimal  approximations.  Us  main  attractions  are 
tlrat  it  is  deterministic  and  efficient. 


7.  Discussion  and  Research  Directions 

Several  issues  concerning  the  framework  for  computing  visible-surface  representations  arc  discussed 
in  tins  section,  and  directions  for  future  research  arc  suggested.  Ihe  discussion  focuses  on 
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Figure  36.  Detected  discontinuities. 


discontinuity  detection,  choosing  constraint  parameters,  handling  rivalries  in  constraints,  grouping 
constraints,  invariance  properties  of  the  surface  reconstruction  model,  and  visible-surface  analysis. 
[Tcr/opoulos,  1984,  Ch.  11]  contains  a  more  extensive  treatment  of  these  and  other  issues,  including 
multircsolution  relative  depth  representations  of  surfaces,  and  the  possibility  of  computing  visible- 
surface  representations  “instantaneously”  by  analog  networks. 

7.1.  On  Discontinuity  Detection 

Some  recent  work  in  image  restoration  is  of  relevance  to  the  problem  of  piecewise  continuous 
surface  reconstruction.  A  piecewise  constant  image  model  employed  by  Blake  [1983]  for  image 
reconstruction  is  interesting  in  diat  it  incorporates  “weak  constraints"  which  can  be  broken  at  a  cost. 
The  resulting  optimization  problem  is  related  to  our  variational  continuity  control  approach,  but 
more  restricted.  Blake  used  an  adaptive  method,  which  he  referred  to  as  “graduated  nonconvexity," 
to  obtain  good  solutions  to  the  nonconvcx  problem.  It  has  not  been  established  however  whether 
this  interesting  method  applies  to  the  sparse  data  case  as  well. 

Gcman  and  Gcman  [1985]  used  Markov  random  field  models  with  associated  Gibbs  distributions 
to  restore  piecewise  constant  images  corrupted  by  additive  Gaussian  noise.  The  restoration  seeks 
a  maximum  a  posteriori  estimate  of  the  original  image,  given  die  degraded  image,  and  includes 
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an  explicit  "line  process"  that  estimates  the  locations  of  step  edges  in  intensity.  This  work  was 
restricted  to  dense  image  data.  The  Gcmans’  approach  was  adopted  with  encouraging  results  to 
surface  reconstruction  from  sparse  depth  constraints  by  Marroquin  [1984],  His  markov  random  field 
model,  while  less  restrictive  than  the  Gcmans'  piecewise  constant  one.  in  fact  models  a  membrane 
spline  whose  smoothness  is  insufficient,  for  computing  visible-surface  representations.  A  line  process 
essentially  equivalent  to  the  Gcmans’  was  incorporated  to  estimate  depth  discontinuities.  The 
numerical  solution  strategy  in  both  of  the  above  studies  was  stochastic  optimization  using  the 
Metropolis  algorithm  and  simulated  annealing  to  optimize  the  nonconvcx  functional  [Kirkpatrick 
el  al.,  1983|.  This  strategy  can  find  optimal  solutions,  but  for  such  large  reconstruction  problems  it 
has  been  observed  to  converge  notoriously  slowly,  based  on  our  experience,  we  believe  that  it  can 
be  accelerated,  perhaps  enough  to  make  it  practical,  through  the  use  of  mulliresolution  processing. 

Obviously,  the  line  processes  used  in  the  above  work  as  well  as  our  own  encoding  of  disconti¬ 
nuity  contour  configurations  is  unpieasingly  heuristic  and  in  need  of  refinement.  The  discontinuity 
map  can  be  augmented  by  nodal  variables  to  encode  the  local  orientations  of  the  curvilinear  ele¬ 
ments  to  a  higher  degree  of  accuracy.  Such  an  encoding  is  employed  by  /.ucker  and  Parent  (1984) 
in  an  optimization  (relaxation  labeling)  approach  to  finding  contours  in  images.  It  appears  that 
ideas  from  their  work  can  also  be  applied  to  finding  surface  discontinuities  within  our  framework. 

A  promising  possibility  is  to  employ  II)  controllcd-continuity  stabilizers  as  formal  models 
of  smoothness  constraints  along  surface  discontinuity  contours  in  the  x-y  plane.  A  functional 
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Figure  38.  Scatlered  dcpih  constraints  consistent  with  sloping  planes  meeting  discontinuous^  (top)  and  the 
smooth  reconstructed  surface  (bottom). 


that  naturally  comes  to  mind  is  the  curvilinear  analog  of  the  thin  plate  surface  under  tension: 
§cPi(s){ftn{s){()2cl()s2)  f  [  1  -  /)<, ( * ) ] ( /) r / <) ) }  ds .  where  s  denotes  arc  length  along  discontinuity 
contours  c  £  C.  Here,  ftb  allows  breaks,  while  /?„  allows  angles  (tanjent  discontinuities)  to  form 
in  the  discontinuity  contours.  Again,  additional  energy  penalties  must  be  associated  with  these 
occurrences.  Given  our  finite  element  representation  of  surfaces,  curvilinear  finite  elements  are 
the  natural  local  representation  for  discontinuity  contours.  The  combined  variational  principle  has 
both  a  surface  component  and  an  analogous  contour  component.  Although  technically  nontrivial, 
a  formulation  of  surface  reconstruction  generalized  along  these  lines  has  very  strong  appeal. 

7.2.  On  Constraints  —  Parameters,  Rivalries,  and  Grouping 

"ITic  constraint  (spring)  parameters  offer  the  flexibility  to  individually  tunc  the  cocrcivcncss  of  each 
constraint  on  the  reconstructed  surface.  In  the  special  case  of  Gaussian  error  distributions,  die 
parameters  should  be  inversely  proportional  to  the  expected  variances  (o,  -  l/Aa?).  It  ought  to  be 
possible  for  the  low-level  visual  processes  to  associate  a  variance  estimate  or  confidence  with  each 
constraint  that  they  provide.  In  general,  however,  it's  not  obvious  how  to  choose  die  constraint 
parameters  optimally. 

The  constant  of  proportionality  A"1  can  also  be  used  to  tune  die  overall  smoothness  of 
the  reconstructed  surface.  Cross  validation  techniques  may  be  used  to  set  A  optimally  (c.g., 
[VVahba  and  Wcndclbcrger,  1980]).  The  basic  criterion  is  to  choose  A  so  as  to  minimize  over  all 
constraints  the  (weighted)  discrepancy  between  each  constraint  and  its  value  as  estimated  from  the 
surface  reconstructed  using  the  remaining  constraints.  Unfortunately,  diis  involves  computationally 
expensive  sequential  algorithms.  Interestingly,  the  continuous  tuning  of  surface  smoothness  is 
analogous  to  the  scale  space  filtering  technique  proposed  by  Wilkin  [1983]  with  the  added  attraction 
that  it  can  be  applied  to  scattered  data. 

Although  the  variational  principle  was  designed  to  account  for  measurement  errors  in  the 
constraints,  the  possibility  of  massive  rivalries  between  constraints  from  different  sources,  such  as 
stcrcopsis  and  analysis  of  motion  processes,  was  disregarded.  Massive  rivalries  are  unnatural  visual 
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Figure  39.  Evolution  of  the  discontinuity  detection  process. 


phenomena  that  can  nevertheless  occur,  especially  under  contrived  conditions,  and  they  often  lead 
to  multistable  percepts  [Attncavc,  1971].  The  framework  can  potentially  accommodate  rivalries  with 
a  mechanism  that  inhibits  individual  or  entire  sets  of  constraints  by  nullifying  selected  constraint 
parameters.  This  mechanism  can  be  activated  by  a  global  arbitrator  which  monitors  the  contents 
of  visible-surface  representations  to  detect  rivalries  may  also  have  access  to  higher  level  knowledge 
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about  the  scene.  The  arbitrator’s  influence  can  account  for  multistability. 

A  particular  type  of  rivalry  arises  from  transparent  surfaces.  For  instance,  a  surface  such  as  a 
dirty  window  in  front  of  a  background  scene  would  lead  to  two  well  defined  populations  of  depth 
(and  orientation)  constraints  over  the  same  visual  angle,  one  from  the  window,  the  other  from  the 
background.  A  transparency  interpretation  can  be  arrived  at  by  an  arbitrator  which  monitors  the 
surface  reconstruction  process  looking  for  high  approximation  error  between  surface  and  constraints 
over  a  significant  area.  Under  the  e  conditions  the  arbitrator  can  trigger  a  constraint  grouping 
process  which  clusters  the  constraints  into  two  populations,  based  on  depth  values,  say.  Multiple 
surfaces  can  then  be  reconstructed  over  the  same  visual  area  for  each  resulting  constraint  population. 
This  scheme  has  been  applied  on  a  transparent  surface  random  dot  stereogram  [Terzopoulos,  1984, 
Ch.  11]. 

7.3.  On  Invariance  Properties  of  the  Surface  Model 


As  a  transformation  from  sparse  constraints  to  dense  surfaces,  the  thin  plate  under  tension  model 
can  be  shown  to  be  invariant  under  (i.e.  commutes  with)  certain  image  plane  transformations 
applied  to  the  constraints;  namely,  translations,  rotations,  and  similarity  transformations.  This 
implies  that  surface  shapes  will  be  preserved  through  rigid  motions  of  the  scene  or  viewpoint 
parallel  to  the  image  plane  or  along  the  view  direction.  These  are  essential  invariance  properties 
for  visible-surface  reconstruction  [Terzopoulos,  1982]. 

Note,  however,  that  the  thin  plate  spline,  characterized  by  the  small  deflection  approximation 
J7<4  +  2 v\y  -i-  v*y  dr  dy  to  the  bending  energy  density  of  a  thin  plate,  is  not  invariant  under 
arbitrary  3D  transformations  of  the  constraints.  Thus,  surface  interpolation  using  this  expression  is 
not  invariant  under  changes  in  the  view  direction,  as  Blake  [1984]  points  out.  He  shows  that  rotating 
the  view  direction  induces  the  ID  analog  of  the  thin  plate  spline  to  “wobble,”  and  he  demonstrates 
that  this  effect  is  most  pronounced  as  the  (continuous)  spline  is  inclined  sharply  with  respect  to  the 
viewer  or  is  forced  to  bend  sharply.  Blake  views  this  as  a  problem  that  should  be  eliminated  by 
employing  the  large  deflection  bending  energy  of  the  thin  plate,  a  convex  combination  of  the  mean 
and  Gaussian  curvatures  of  the  surface  v(x,y),  which  is  view  direction  invariant. 

Although  £/>r(t;)  can  also  be  made  view  direction  invariant  by  employing  the  large  deflection 
counterparts  for  the  thin  plate  and  membrane  bending  energies,  this  approach  has  a  serious 
technical  drawback  [  Terzopoulos,  1984]:  The  large  deflection  formulas  lead  to  an  extremely  difficult 
nonlinear  problem  (c.g„  the  large  deflection  equations  for  the  thin  plate  arc  two  coupled  nonlinear 
fourth-order  partial  differential  equations  known  as  Von  Karmann’s  equations  [Szilard,  1974]). 

Fortunately,  the  surface  reconstruction  model,  as  it  stands,  is  not  hampered  by  die  lack  of 
view  direction  invariance  because  the  available  constraints  arc  usually  sufficiently  dense  in  practice 
to  tightly  determine  surface  shape;  as  the  view  direction  is  varied,  the  reconstructed  surface  would 
vary  negligibly  (note  that  Blake’s  experiments  reveal  a  significant  wobble  effect  just  in  the  case 
of  extremely  sparse  constraints),  ITirdicnnoro.  the  explicit  introduction  of  depth  and  orientation 
discontinuities  alleviates  much  of  the  wobble  precisely  at  those  places  where  Blake’s  experiments 
show  it  to  be  most  pronounced  on  a  globally  continuous  surface.  An  interesting  psychophysical 
experiment  would  be  to  determine  whether  there  might  be  some  slight  variance  in  the  surfaces 
perceived  by  humans  viewing  sparse  random  dot  stereograms  while  llie  dots  undergo  simulated 
rigid  31)  transformations  and,  if  so,  whether  die  variations  arc  consistent  with  die  reconstruction 
model  (J.  Mayhcw,  personal  communication). 
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7.4.  On  Visible-Surface  Analysis 

The  visible-surface  representation  is  an  intermediate  and  volatile  description  of  the  30  surfaces 
in  scenes.  It  drives  ensuing  processes  which  generate  stable  higher-level  representations  of  shape 
that  are  better  tuned  to  object  recognition.  The  processing  begins  with  visible-surface  analysis 
whose  goal  is  to  abstract  from  the  numeric,  viewer-centered  representation  a  rich  set  of  more 
symbolic,  object-centered  features  that  are  stable  through  viewpoint  changes.  The  extraction  of 
geometric  surface  features  is  facilitated  by  the  dense  shape  information  provided  by  visible-surface 
representations. 

A  promising  approach  to  visible-surface  analysis  is  to  apply  concepts  from  differential  geometry 
[do  Carmo,  1976).  For  instance,  a  surface’s  intrinsic  geometry  (including  Gaussian  curvature, 
geodesics,  etc.)  is  determined  completely  by  the  first  fundamental  form,  which  defines  arc  length 
over  the  surface.  It’s  extrinsic  geometry  (including  normal  curvature,  principal  curvatures,  etc.)  are 
determined  by  the  second  fundamental  form,  which  describes  the  deviation  of  the  surface  from  the 
local  tangent  plane.  The  fundamental  theorem  of  the  local  theory  of  surfaces  (usually  attributed 
to  Bonnet)  states  that  the  analytic  study  of  surface  properties  consists  of  the  study  of  lire  two 
fundamental  forms;  i.e.,  the  six  fundamental  tensor  coefficients  (which  arc  not  all  independent)  as 
functions  of  the  two  independent  parameters  of  the  surface.  The  fundamental  forms  arc  invariant 
under  changes  in  the  parameterization,  and  together  they  determine  surface  shape  up  to  rigid  body 
transformations.  These  properties  make  them  ideal  foundations  for  object  centered  symbolic  surface 
representations. 

The  visible-surface  representation  makes  it  possible  to  estimate  the  first  and  second  funda¬ 
mental  forms  on  a  point-by-point  basis  over  the  entire  visible  surface.  The  finite  element  shape 
representation  reduces  the  computation  of  crucial  local  surface  features  such  as  die  Gaussian  curva¬ 
ture,  principal  curvatures,  and  principal  directions  to  the  evaluation  of  simple  algebraic  expressions 
of  neighboring  nodal  variables  (see  [Tcrzopoulos,  1984,  Ch.  11]  for  derivations).  It  is  then  a  simple 
step  to  determine  the  elliptic,  hyperbolic,  parabolic,  umbilic,  and  planar  points,  as  well  as  geodesics, 
asymptotes,  and  lines  of  curvature. 

For  example.  Fig.  39  shows  the  reconstructed  surface  of  a  lightbulb.  Fig.  40  shows  the  Gaussian 
curvature  K(x,y)  computed  for  the  reconstructed  lightbulb  surface  of  Fig.  20.  The  elliptic  points 
(K  >  0)  are  shown  in  white,  the  hyperbolic  points  ( K  <  0)  arc  shown  in  black,  and  the  parabolic 
(K  =  0)  points  separate  the  two  regions.  Note  the  alternation  in  the  sign  of  curvature  at  die 
screw  mount.  Fig.  41  plots  the  computed  field  of  principal  directions  for  die  lightbulb  at  the 
two  coarsest  scales.  These  demonstrations  illustrate  die  feasibility  of  reliably  computing  from  these 
representations  higher-order  intrinsic  and  extrinsic  properties  of  surface  shape.  The  reliability  can 
be  attributed  to  the  regularizing  properties  of  the  thin  plate  surface  under  tension  which  overcomes 
the  potentially  detrimental  effects  of  noise  in  the  data,  while  preserving  discontinuities.  For  further 
analysis  of  the  kinds  of  features  that  can  be  computed  from  dense,  numeric,  representations  of 
surfaces  see,  e.g.,  [Brady  el  al„  1985]  or  [Mcdioni  and  Ncvatia,  1984]. 


8.  Conclusion 

Constraints  on  surface  shape,  contributed  by  multiple  low-level  visual  processes,  can  be  computed 
reliably  at  multiple  resolutions,  but  only  at  scattered  locations  in  the  field  of  view.  Subsequent 
visual  processing  can  be  facilitated  substantially  if  the  scattered  constraints  are  transformed  into 
visible-surface  representations  that  make  surface  shape  explicit  everywhere.  To  accomplish  this 
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Figure  40.  Elliptic  (white)  and  hyperbolic  (black)  points  of  the  reconstructed  lightbulb  at  four  scales. 


effectively,  information  must  be  integrated  over  multiple  visual  modalities  and  fused  across  multiple 
scales  of  resolution. 

In  this  paper,  we  have  developed  a  computational  theory  of  visible-surface  representations. 
Within  a  unified  computational  framework,  formal  solutions  were  offered  to  fundamental  problems 
of  reconstructing  visible  surfaces:  (i)  integrating  constraints  on  the  depth  and  orientation  of  surfaces 
across  various  modalities  and  scales,  (ii)  interpolating  surface  shape  information  into  (piecewise) 
smooth  surfaces,  (iii)  discovering  discontinuities  in  surface  depth  and  orientation  and  enabling  them 
to  restrict  interpolation,  and  (iv)  efficiently  maintaining  consistency  in  distributed,  multircsolution 
visible-surface  representations. 

A  visible-surface  reconstruction  algorithm  implements  the  framework.  Extensive  testing  has 
shown  it  to  be  viable.  'Ihe  algorithm  coordinates  cooperative  processes  within  a  multircsolution 
hierarchy  of  surface  representations  to  dramatically  increase  computational  efficiency.  It  is  well 
suited  to  implementation  on  massively  parallel  networks  of  simple,  locally  interconnected  processors. 
Such  computational  networks  arc  suggestive  of  biological  mechanisms  and  are  also  well  suited  to 
VLSI  technology. 
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Figure  41.  Principal  directions  for  the  reconstructed  lighlbuib  at  the  two  coarsest  scales.  The  directions  of 
greatest  curvature  are  shown  on  the  left.  Those  of  least  curvature  are  shown  on  the  right 
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